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Abstract 

A  numerical  method  with  the  approximate  global  convergence  property  is  developed 
for  a  3-D  Coefficient  Inverse  Problem  for  a  hyperbolic  PDE  with  the  backscattering 
data.  An  important  part  of  this  technique  is  the  quasi-reversibility  method.  Approxi¬ 
mate  global  convergence  theorem  is  proven.  Results  of  two  numerical  experiments  are 
presented. 
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1  Introduction 

The  goal  of  this  publication  is  to  extend  the  numerical  method  of  [4,  5,  6,  7,  8,  27]  to  the  case 
of  backscattering  data  resulting  from  a  single  measurement.  “Single  measurement”  means 
that  the  time  dependent  backscattering  data  for  a  Coefficient  Inverse  Problem  (CIP)  for  a 
hyperbolic  PDE  are  generated  by  a  single  location  of  the  point  source.  We  develop  a  new 
theory  and  present  numerical  results  confirming  this  theory.  Our  target  application  is  in 
imaging  of  plastic  antipersonnel  land  mines,  and  we  model  this  application  in  our  numerical 
testing.  In  military  applications,  due  to  various  dangers  on  the  battlefield,  it  is  desirable  to 
minimize  both  the  number  of  measurements  and  the  observation  angle.  In  the  current  paper 
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we  model  the  most  suitable  arrangement  for  this  case,  which  is  to  use  a  single  position  of 
the  point  source  and  to  measure  only  the  backscattering  signal. 

It  will  be  shown  in  the  forthcoming  publication  [31]  that  a  close  analog  of  the  numerical 
method  of  this  paper  provides  accurate  results  for  the  most  challenging  case  of  blind  backscat¬ 
tering  experimental  data  collected  from  targets  mimicking  antipersonnel  land  mines.  These 
data  were  collected  in  the  field  by  a  new  forward  looking  radar,  see  [36]  for  a  description  of 
this  radar. 

This  is  a  significantly  revised  version  of  the  paper  [30].  Main  analytical  results  of  the 
current  paper:  Theorems  3.1,  6.1,  7.1,  Lemma  7.1  and  topics  of  subsections  6.2  and  6.3  are 
not  parts  of  [30].  Lemmata  5.1,  5.2  were  formulated  but  not  proven  in  [30].  Lemma  5.3  was 
proven  in  [30].  Theorem  5.1  was  proven  in  [30]  for  the  2-D  case.  Since  the  proof  in  3-D 
is  similar,  we  omit  it  here.  Numerical  examples  of  this  paper  as  well  as  Figures  8. 1,8. 2  are 
different  from  ones  of  [30]. 

1.1  Preliminaries 

It  is  well  known  that  the  goal  of  the  construction  of  reliable  numerical  methods  for  CIPs 
faces  quite  substantial  challenges.  They  are  mainly  caused  by  two  factors  combined:  non¬ 
linearity  and  ill-posedness  of  CIPs.  Least  squares  regularizing  Tikhonov  functionals  suffer 
from  the  commonly  known  phenomenon  of  multiple  local  minima  and  ravines.  Therefore,  a 
gradient-like  method  of  the  minimization  of  this  functional  can  stop  at  any  point  of  either 
a  local  minimum  or  a  ravine.  However,  this  point  can  be  located  quite  far  from  the  true 
solution.  This  results  in  the  local  convergence  of  conventional  numerical  methods  for  CIPs, 
such  as,  e.g.  various  versions  of  the  Newton  and  gradient  methods.  Convergence  of  these 
algorithms  is  guaranteed  only  if  the  starting  point  of  iterations  is  located  in  a  sufficiently 
small  neighborhood  of  the  exact  solution.  At  the  same  time,  in  many  applications  a  good 
approximation  for  the  exact  solution  is  rarely  known  in  advance. 

Therefore  the  central  question  in  a  computational  treatment  of  a  CIP  is  to  develop  such 
a  numerical  method,  which  would  provide  at  least  one  point  in  a  sufficiently  small  neighbor¬ 
hood  of  the  exact  solution  without  any  a  priori  knowledge  of  that  neighborhood.  Further¬ 
more,  the  latter  property  should  be  rigorously  guaranteed,  at  least  within  the  framework  of  a 
certain  reasonable  approximate  mathematical  model.  In  addition,  numerical  studies  should 
confirm  this  property. 

This  question  was  addressed  in  a  series  of  recent  publications  [4,  5,  6,  7,  8,  27]  for  CIPs 
for  a  hyperbolic  PDE  in  the  case  when  the  data  resulting  from  a  single  measurement  are 
given  at  the  entire  boundary,  i.e.  the  case  of  complete  data  collection.  The  technique  of 
neither  these  publications  nor  the  current  one  does  not  use  least  squares  functionals.  Only 
the  structure  of  the  underlying  PDE  operator  is  used.  Also,  it  does  not  use  a  knowledge  of 
the  background  values  of  the  unknown  coefficient. 

Because  of  a  quite  significant  difficulty  of  addressing  the  above  question,  it  is  unlikely 
that  this  can  be  done  without  some  approximations.  Hence,  a  reasonable  approximate  math¬ 
ematical  model  was  proposed  on  pages  102-104  of  [5]  and  on  page  11  of  [27].  Developing  this 
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concept  further,  we  introduce  in  subsection  6.2  a  new  concept  of  the  so-called  approximate 
global  convergence  property.  Based  on  this  concept,  we  specify  in  subsection  6.3  our  new 
approximate  mathematical  model. 

As  soon  as  the  desired  first  good  approximation  is  obtained,  a  locally  convergent  numer¬ 
ical  method  can  be  applied  to  refine  it,  see,  e.g.  books  [2,  24]  for  Newton-like  methods  for 
ill-posed  problems.  Thus,  a  two-stage  numerical  procedure  was  developed  in  [5,  6,  7,  8,  30], 
and  it  also  takes  place  here.  On  the  first  stage  a  good  first  approximation  is  provided  by  the 
above  technique.  On  the  second  stage  this  approximation  is  refined  via  a  locally  convergent 
method.  This  method  takes  the  first  stage  solution  as  the  starting  point  for  iterations. 

In  addition  to  the  synthetic  data  of  numerical  studies  of  [4,  5,  6,  7],  an  experimental 
verification  of  the  first  stage  was  carried  out  in  [27]  for  the  case  of  blind  experimental  data 
collected  in  the  transmitted  scattering  mode.  A  very  accurate  reconstruction  of  both  re¬ 
fractive  indices  and  locations  of  dielectric  inclusions  was  demonstrated  in  [27].  Next,  the 
two-stage  procedure  was  verified  on  the  same  experimental  data  in  [8].  A  very  good  re¬ 
construction  accuracy  of  locations,  shapes  and  refractive  indices  of  dielectric  inclusions  was 
observed  in  [8]. 

1.2  What  is  new  in  this  paper 

Since  we  have  over-determined  boundary  conditions  on  the  backscattering  part  of  the  bound¬ 
ary  of  the  domain  of  interest,  we  use  the  Quasi- Reversibility  Method  (QRM)  [28,  33].  The 
QRM  is  well  suited  for  this  kind  of  problems,  since  it  finds  the  “least  squares”  solution. 
The  QRM  was  not  in  [4,  5,  6,  7,  8,  27]  since  the  complete  data  collection  was  used  in  these 
references.  We  refer  to,  e.g.  [14,  15,  16,  18,  25]  for  some  publications  on  the  QRM  for  linear 
ill-posed  problems.  The  convergence  analysis  of  the  QRM  is  based  on  Carleman  estimates. 

The  goal  of  the  analytical  part  of  this  paper  is  to  establish  the  above  mentioned  approx¬ 
imate  global  convergence  property  for  our  algorithm.  While  the  QRM  is  well  known  for 
linear  ill-posed  problems,  its  application  to  a  nonlinear  problem,  such  as  our  CIP,  is  new. 
The  only  exception  in  this  regard  is  the  recently  published  work  [29]  of  the  first  and  third 
authors,  where  the  technique  of  [4,  5]  was  extended  to  the  1-D  case  with  the  backscattering 
data  and  the  QRM  was  applied.  However,  the  Carleman  estimate  in  1-D  is  simpler  than  the 
one  in  n— D  (n  =  2,  3)  since  in  1-D  it  enables  to  estimate  the  H2— norm  of  the  solution  in 
the  whole  interval  of  interest.  Unlike  this,  the  Carleman  estimate  in  the  n— D  (n  =  2,  3)  case 
enables  one  to  obtain  the  Holder  stability  estimate  of  the  H2— norm  of  the  QRM  solution 
only  in  a  subdomain  of  the  domain  of  interest.  This  causes  the  first  main  difficulty  in  the 
convergence  analysis  compared  with  the  linear  case. 

Let  a  G  (0, 1)  be  a  small  regularization  parameter  of  the  QRM.  The  desired  estimate 
of  the  H2 —norm  of  the  solution  in  a  subdomain  inevitably  contains  the  large  multiplier 
a~r,r  =  const.  >  0.  This  multiplier  was  not  a  problem  in  the  linear  case  since  the  QRM 
is  applied  only  once  in  this  case.  Now,  however,  because  of  the  nonlinearity  of  our  CIP, 
we  need  to  arrange  an  iterative  process.  Thus,  we  need  to  “suppress”  this  multiplier  on 
each  iteration  to  ensure  stability.  This  causes  the  second  main  difficulty,  compared  with  the 
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linear  case.  The  third  main  difficulty  here,  compared  with  the  linear  case,  is  that  along  with 
the  if2— norm  of  the  QRM  solution  in  the  subdomain,  we  need  to  iteratively  estimate  its 
H 5—  norm  in  the  entire  domain  (because  of  the  regularization  term),  and  the  latter  estimate 
contains  a  large  multiplier  a_1//2  which  again  needs  to  be  “suppressed”  on  each  iteration  to 
stabilize  the  process.  Four  more  new  analytical  elements  are: 

1.  To  perform  our  convergence  analysis,  we  establish  in  Theorem  3.1  (section  3)  some 
new  properties  of  the  solution  of  an  elliptic  PDE  in  M3,  including  estimates  of  this  solution 
both  from  the  above  and  from  the  below.  This  PDE  is  obtained  via  the  Laplace  transform 
of  the  originating  hyperbolic  PDE.  The  difficulty  here  is  that  while  the  classical  theory  of 
elliptic  PDEs  is  developed  only  for  bounded  domains  [21],  we  need  to  work  in  M3. 

2.  We  obtain  accuracy  estimates  for  the  so-called  “tail”  functions  on  each  iteration.  Such 
estimates  were  not  obtained  in  above  cited  previous  publications  about  this  method.  These 
estimates  lead  to  a  more  convenient  (than  previously)  formulation  of  the  approximate  global 
convergence  Theorem  7.1. 

3.  We  use  a  new  Carleman  Weight  Function  (CWF)  in  our  Carleman  estimate  for  the 
Laplace  operator.  This  CWF  is  well  suitable  for  a  boundary  condition,  which  is  used  here 
for  a  better  stability  of  the  QRM. 

4.  In  the  above  cited  previous  publications  about  this  method  functions  qn^  (section 
4)  were  strong  solutions  of  the  Dirichlct  boundary  value  problems  for  corresponding  elliptic 
PDEs,  and  Schauder  theorem  was  applied  to  estimate  them.  Now,  however,  each  function 
qn.k,  is  a  weak  solution  of  a  certain  integral  identity  since  it  is  the  minimizer  of  the  Tikhonov 
functional  corresponding  to  the  QRM.  This  causes  additional  difficulties  in  the  proof  of 
Theorem  7.1. 

While  we  work  with  the  data  resulting  from  a  single  measurement,  we  refer  to  works 
[1,  9,  10,  11,  23,  37]  and  references  cited  there  for  some  non-locally  convergent  algorithms 
for  CIPs  with  the  data  resulting  from  multiple  measurements.  In  particular,  publications 
[9,  10,  23]  consider  CIPs  for  hyperbolic  PDEs.  In  section  2  we  pose  the  inverse  problem.  In 
section  3  we  establish  some  properties  of  the  Laplace  transform  of  the  solution  of  the  forward 
problem.  The  algorithm  is  described  in  section  4.  Sections  5-7  are  about  the  convergence 
analysis.  Numerical  results  are  presented  in  section  8. 


2  Posing  the  Inverse  Problem 

Below  x  =  ( x ,  y,  z )  €  R3  denotes  an  arbitrary  vector  as  well  as  the  first  component  of  this 
vector.  The  context  does  not  allows  an  ambiguity.  As  the  forward  problem,  we  consider  the 
Cauchy  problem  for  the  following  hyperbolic  PDE 

c(x)utt  =  Am  in  R3  x  (0,  oo) ,  (1) 

u  (x,  0)  =  0,  ut  (x,  0)  =  8  (x  —  x0) .  (2) 


Equation  (1)  governs,  e.g.  propagation  of  acoustic  and  electromagnetic  waves.  In  the  acous¬ 
tical  case  1  /  \fcjx)  is  the  sound  speed.  In  the  2-D  case  of  EM  waves  propagation  in  a 
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non-magnetic  medium  c(x)  :=  er  (x) ,  where  er  (. x )  is  the  spatially  distributed  dielectric  con¬ 
stant,  see,  e.g.  [17]  for  the  derivation  of  (1)  from  the  Maxwell  equations  in  the  2-D  case. 
Even  though  equation  (1)  cannot  be  derived  from  the  Maxwell  equations  in  3-D  for  the  case 
er  (. x )  7^  const.,  still  it  was  successfully  used  in  [8,  27,  31]  to  model  the  propagation  of  EM 
waves  for  experimental  data,  which  are  in  3-D. 

Let  fl  C  IR3  be  a  convex  bounded  domain  with  the  piecewise  smooth  boundary  <9fh  As 
it  is  always  the  case  in  CIPs,  we  need  to  assume  a  certain  over-smoothness  of  the  unknown 
coefficient.  Let  d  >  1  be  a  certain  given  number,  which  is  the  upper  bound  for  our  unknown 
coefficient  c.  Note  that  some  a  priori  information  about  the  unknown  coefficient  should  be 
given  in  accordance  with  the  Tikhonov  concept  for  Ill-Posed  Problems  [40].  However,  since 
we  do  not  assume  a  smallness  of  the  number  d  —  d  —  1,  then  a  smallness  assumption  is  not 
imposed  on  the  function  c.  So,  we  assume  that 

c  (x)  E  [1,  d] ,  c  (x)  —  1  for  x  E  M3\h2,  (3) 

c(x)  E  C4(M3).  (4) 

The  assumption  that  c  (x)  is,  in  general,  less  outside  of  the  domain  of  interest  D  than  inside 
of  it  is  because  in  our  target  application  to  imaging  of  antipersonnel  plastic  land  mines 
the  dielectric  constant  of  mine-like  targets  is  greater  than  the  one  in  the  background,  see 
subsection  8.2.  To  simplify  the  presentation  and  also  because  of  our  target  application,  we 
now  specify  the  domain  fl  C  M3  as  follows.  Let  P  >  0  be  a  constant.  Below 


Q 

=  {{x,y,z) 

:  -P  < 

x,y  <  P,z  E  (0,  2P)}  ,  dn  =  U3=1Tj, 

(5) 

r\ 

=  {(x,y,z) 

:  -P  < 

x,y  <  P,z  =  0}  , 

(6) 

r2 

=  {(x,y,z) 

:x,y  = 

±P,  z  E  (0,  2P)}  ,  r3  =  n  n  {z  =  2 P}  . 

(7) 

Inverse  Problem.  Suppose  that  the  coefficient  c  (x)  in  equation  (1)  satisfies  conditions 
(3),  (4)  and  is  unknown  in  the  domain  fb  Determine  c(x)  for  x  E  D,  assuming  that  the 
following  functions  Tp0  (x,  t )  and  Tp1  (x,  t )  are  known  for  a  single  source  position  xo  f  D 

u  (x,  t)  |n  =  Od  t) ,  uz  (x,  t )  |n  =  (x,  t)  ,t  e  (0,  oo) .  (8) 

In  experiments  usually  only  the  function  TpQ  (x,  t)  is  measured.  One  can  approximately 
assume  that  this  function  is  known  at  the  entire  plane  {z  =  0} .  Next,  since  by  (3)  and  (5) 
the  coefficient  c(x)  —  1  for  z  <  0,  then  solving  the  forward  problem  (1),  (2)  in  the  half  space 
{z  <  0}  with  the  boundary  condition  u  (x,  t )  \z=o—  {xi  £)>  one  can  uniquely  determine  the 
function  u  (x,  y ,  z,  t )  for  z  <  0,  t  >  0,  which  gives  Tp1  (x,  t ). 

The  next  reasonable  question  to  address  is  about  the  infinite  rather  than  a  finite  time 
interval  at  which  functions  Tp0  (x,  t )  ,Tpx  (x,  t)  are  given  in  (8).  In  our  numerical  procedure  we 
use  the  Laplace  transform  with  respect  to  t.  Since  the  kernel  e~st,  s  >  0  of  this  transform 
decays  rapidly  when  t  — >  oo,  then  the  integral  over  a  finite  time  interval  is  approximately 
the  same  as  the  one  over  the  infinite  interval.  This  justifies  our  assumption  t  E  (0,  oo)  in 
(8).  In  addition,  when  working  with  experimental  data,  which  were  measured  on  a  finite 
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time  interval  [8,  27,  31],  we  have  observed  that  the  assumption  t  G  (0,  oo)  does  not  affect 
the  quality  of  solutions. 

Remark  2.1.  The  question  of  the  uniqueness  of  this  CIP  is  a  well  known  long  standing 
problem.  Currently  it  is  addressed  positively  via  the  method  of  Carleman  estimates  only 
in  the  case  when  5  (x  —  x0)  in  (2)  is  replaced  with  a  function  f  (x)  ^  0  in  12  [26,  28].  Still, 
because  of  applications,  it  makes  sense  to  develop  numerical  methods  for  this  CIP.  Thus,  we 
assume  everywhere  below  that  uniqueness  of  our  Inverse  Problem  holds,  also  see  Lemma  6.1. 


3  Laplace  Transform 

Below  Ck+1  are  Holder  spaces,  where  k  >  0  is  an  integer  and  7  G  (0, 1).  Let  the  function 
c  ( x )  be  such  that 

c  G  C1  (R3)  and  conditions  (3)  hold.  (9) 

Consider  the  Laplace  transform  w  of  the  function  u.  We  are  not  concerned  with  inverting 
this  transform,  since  we  use  it  only  to  approximate  the  unknown  coefficient  c  ( x ) .  Thus, 

OO 

for  s  >  s  —  const.  >  0,  (10) 

0 

where  s  >  0  is  a  certain  number.  In  our  numerical  studies  we  choose  s  experimentally.  We 
call  the  parameter  s  pseudo  frequency.  The  function  w  is  satisfies  the  following  conditions 

A w  —  s2c(x)w  =  —  5  (x  —  x0)  ,x  G  R3,  Vs  >  s,  (11) 

lim  w(x,s)  =  0,Vs>s.  (12) 

|jc|— kx> 


The  condition  (12)  was  established  in  subsection  2.2  of  [6]  for  s  >  s°,  where  the  number 
s°  =  s°  (j|c||C7(R3)j  >  0  is  sufficiently  large.  It  was  also  proven  in  [6]  that  for  s  >  s°  the 
function  w  (. x ,  s)  has  the  form 


w 

w0  (x,  s ) 


w0  +  w,  where  w  G  C2+7  (R3)  , 
exp  (— s  \x  —  x0|) 

47T  \x  —  Xq\ 


(13) 

(14) 


The  function  wq  solves  the  problem  (11),  (12)  for  the  case  c(x)  =  1.  Theorem  6.17  of  [21] 
about  the  smoothness  of  solutions  of  elliptic  PDEs  gives  w  G  Ck+2+1  (R3) ,  if  c  G  Ck+1  (R3) . 

In  our  algorithm  as  well  as  in  the  course  of  the  proof  of  Theorem  7.1  we  will  construct 
solutions  w(x,s)  :=  w(x,s;c )  of  the  problem  (11),  (12)  with  a  certain  s  >  1  for  different 
functions  c(x)  satisfying  (9).  These  functions  c(x)  will  be  calculated  iteratively.  We  will 
not  estimate  C 7  (R3)  —  norms  of  those  functions  c(x).  Hence,  if  constructing  the  sequence 
of  functions  w(x,s‘,c )  via  solving  the  forward  problem  (1),  (2)  and  then  calculating  the 
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integral  (10)  for  s  :=  s,  then  each  of  those  functions  c  (x)  might  require  its  own  value 
s°  :=  s°  (j|c||C7(R3))  ,  and  there  is  no  guarantee  that  s°  (j|c||C7(1R3) j  <  s.  However,  since  we 
will  work  only  with  the  interval  s  G  [s,  s]  for  fixed  values  of  parameters  s,  s,  then  we  need  to 
ensure  that  for  each  such  function  c  and  for  each  value  s  >  0  there  exists  unique  solution  of 
the  form  (13),  (14)  of  the  problem  (11),  (12).  Thus,  we  now  prove 

Theorem  3.1.  Let  x0  ^  O,  the  function  c(x)  satisfies  condition  (9)  and  also  c  G 
Ck+i  (]g>3)  _  Tfien  jor  any  s  >  o  there  exists  unique  solution  of  the  problem  (11),  (12)  of  the 
form  (13),  (If)  with  w  G  Ck+2+1  (R3) .  Denote 


( x ,  s)  = 


exp 


iy/d\ 


x  —  x0 


An  \x  —  x0\ 

the  solution  of  the  problem  (11),  (12)  for  c(x)  =  d.  Then 

(x,  s)  <  w  (x,  s )  <  wq  (x,  s) ,  Vx  7^  Xq- 
Proof.  Consider  the  following  parabolic  Cauchy  problem  for  (x,  t)  e 


(15) 


(16) 


X  (0,  oo) 


c  (x)  vt  =  An,  v  (x,  0 )  —  S  (x  —  x0)  ■ 


(17) 


Let  v0(x,t )  be  the  solution  of  the  problem  (17)  with  c  =  1.  Also,  consider  the  function 
v  (x,  t) , 


v0  (x, t ) 


x  -  x0|2  \ 

At  J  ’ 


(18) 


t 

(■ v  —  Vo)  (x,  r )  dr. 
o 

Denote  b(x)  —  c(x)  —  1.  We  obtain  from  (17)  and  (19) 

Av  —  c  (x)  vt  =  b  (x)  v0,  v  (. x ,  0)  =  0,  (x,  t)  G  IR3  x  (0,  oo) . 


(19) 


(20) 


Since  x0  ^  and  by  (3)  and  (9)  b  (x)  —  0  outside  of  then  it  follows  from  (18)  that 
the  right  hand  of  (20)  does  not  have  a  singularity  in  R3  x  [0,oo).  Let  T,R  >  0  be  two 
arbitrary  numbers  and  BR  (T)  =  {|x|  <  R}  x  (0,  T ) .  By  (3),  (9)  and  (18)  b  (. x )  v0  (. x ,  t)  >  0 
for  (x,  t )  G  R3  x  (0,  oo) .  Hence,  applying  to  (20)  the  maximum  principle  of  Theorem  1  of 
Chapter  2  of  [20],  we  obtain  ma x-^R^v  (x,t)  <  0.  Since  numbers  R,  T  >  0  are  arbitrary, 
then 

v  (x,  t)  <  0  in  R3  x  [0,  oo) .  (21) 

On  the  other  hand,  Theorem  11  of  Chapter  2  of  [20]  ensures  that  the  fundamental  solution 
of  the  parabolic  equation  is  positive  for  t  >  0.  Hence,  (19)  and  (21)  imply  that 


v  (x,  t)  dr  <  /  Vo  (x,  t )  dr  and  v  (x,  t )  >  0. 


(22) 


Next,  consider  the  following  version  of  the  Laplace  transform 


C2v  —  I  v(x,t)e  stdt. 
o 

By  one  of  well  known  properties  of  the  Laplace  transform 

t 


(23) 


C, 


f  (t)  dr  =  —C2f 


(24) 


for  any  appropriate  function  /.  Formula  (28)  of  section  4.5  of  Tables  [3]  gives  C2Vq  =  wo- 
Fix  a  number  s  >  0.  Hence,  (22)-(24)  and  Fubini  theorem  lead  to 
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C2  \  I  v  ( X ,  r)  dr  1  =  —C2v  <  —  C2  {v0)  =  —w0  (x,  s) . 

IS  s  s 

,0 

Hence,  the  integral  (23)  converges  absolutely.  Next,  by  (20)  for  any  A  >  0 

A  A  A 

A  J  v  (x,t)  e~s2tdt  —  I  A v  (x,t)  e~s2tdt  =  j  [cvt  +  (c  —  1)  v0]  e~s2tdt. 
00  0 

Setting  here  A  — >  00  and  using  that  by  (19)  cvt  +  (c  —  1)  vQ  =  cv  —  v0,  we  obtain 

A  A 

lim  A  v  (x,  t )  e~s2tdt  =  lim  /  Av  (x,  t )  e~s2tdt  =  cC2v  —  C2Vq. 

A^oo  /  A— xx)  / 


(25) 


(26) 


Hence,  it  follows  from  (26)  that  A C2  (v)  and  C2  (An)  exist  and  A C-2  (v)  =  C2  (Av) .  Further¬ 
more,  by  (19)  and  (24)-(26)  A C2  (v)  =  s^2A  (C2v  —  C2v 0)  =  cC2v  —  C2v 0.  Hence,  denoting 
w  :=  C2  (v)  and  using  C2v0  =  w0  as  well  as  Aw0  —  s2w0  =  —  S  (x  —  Xo) ,  we  obtain  that  the 
function  w  satishes  equation  (11). 

We  now  prove  (12)  and  also  that 


w  =  Wo  +  w,  where  w  G  Ck+2+1  (M3)  . 
Since  cvt  =  vt  +  bvt,  then  using  (19)  and  (20),  we  obtain 

vt  =  Av  —  b  (x)  v ,  v  (x,  0)  =  0. 


(27) 


(28) 


Since  b  e  Ck+1  (M3)  and  b(x)  —  0  near  x0,  then  at  least  v  G  C2+7,1+7//2  (R3  x  [0,  T])  ,VT  >  0 
(formula  (13.2)  of  Chapter  4  of  [32]).  Hence,  by  (27) 


V  (x,  t)  =  -  j  J  v0  (x-€,t-r)b  (0  v  (f ,  r)  d^dr. 
0  n 


(29) 
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By  (19),  (23)  and  (25)  C^v  =  s  2  {C-^v  —  wq)  —  s  2  (w  —  wq ) .  Hence,  applying  the  Laplace 
transform  Co  to  both  sides  of  (29)  and  using  the  convolution  theorem,  we  obtain 


w  ( X ,  s)  =  w0  (x,  s)  -  s2  J  w0  (x-€)b  (0  w  (f ,  s)  d£. 

Cl 


(30) 


Iterating  this  integral  equation  once,  we  obtain 


w  (x,  s )  =  w0  (x,  s)  -  s2  J  w0  (x  -  0  wo  (0  b  (£)  d£ 

n 


+s4  J  w0(x-£)b  (0  J  w0(€-  rj)  b  (rj)  w  (77,  s)  drj.  (31) 

O  Cl 

Since  b  G  C7  (O)  ,  then  the  well  known  results  for  integral  equations  with  kernels  like 
Wq  (x  —  £)  imply  that 


By  (14)  and  (30)  functions  w  (x,  s) ,  Wo  (x,  s)  ,w(x,  s)  =  (w  —  Wo)  (x,  s)  satisfy  condition 
(12).  Furthermore,  it  follows  from  (31)-(33)  that  (27)  holds  with  k  —  0.  Next,  the  above 
mentioned  Theorem  6.17  of  [21]  ensures  that  (27)  is  true  with  k  >  0.  Thus,  we  have  proven 
existence.  Uniqueness  is  proven  in  subsection  2.2  of  [6]. 

Finally  we  prove  (16).  Since  w  =  C2  (v)  and  v  >  0,  b  >  0,  then  the  right  inequality  (16) 
follows  from  (30).  Consider  the  function  w ^  (x,  s)  =  w  (. x ,  s)  —  (. x ,  s)  .  Then  (11),  (12) 

and  (15)  imply  that 

Aw^  —  s2cww  =  s 2  (c  (x)  —  d )  w^d\  lim  w ^  (x,  s )  =  0.  (34) 

|jc|— >00 


By  (13)-(15) 


w^  (x,  s ) 


=  exp 


y/d  —  1 


F  -  ^0 


tfU)  (a;,  5) 

Hence,  there  exists  a  sufficiently  small  number  £  >  0  such  that 

w(cb  (x,  s)  >  0  for  x  G  {|x  —  x0|  <  e}  . 


(1  -  O  (|x  -  x0|))  >  0,  x  ->  a;0,  x  ^  x0. 


(35) 
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For  R  >  0  consider  the  domain  Brj6  =  {|x|  <  R,  \x  —  xq\  >  e} .  Assuming  that  Bri£  ^  0, 
which  is  true  for  sufficiently  large  R,  we  obtain  w(d>  e  C2+1  (B^e)  and  s 2  (c  (x)  —  d )  <  0 

in  Bre.  Hence,  applying  the  maximum  principle  to  equation  (34),  we  obtain  min-^;  wd  > 
min dBReWd-  Setting  R  — >  oo  and  using  the  second  condition  (34)  as  well  as  (35),  we  obtain 
min|x_xo| >£w(d)  >  min|x_xo|=e  wd  >  0.  Thus,  w  (x,s)  >  w (d)  (x,  s)  for  x  x0.  □ 

To  justify  our  new  approximate  mathematical  model  of  subsection  6.3,  we  need  the 
asymptotic  behavior  of  the  function  w(x,s )  at  s  — >  oo,  see  discussion  in  subsection  8.1 
about  the  verification  of  conditions  of  Lemma  3.1. 

Lemma  3.1  [4],  Assume  that  conditions  (3)  and  (4)  hold.  Let  u(x,t )  be  the  solution  of 
the  problem  (1),  (2)  and  for  sufficiently  large  values  of  s  >  0  let  w(x,s )  =  C\u.  Assume 
that  geodesic  lines,  generated  by  the  eikonal  equation  corresponding  to  the  function  c  (x)  are 
regular,  i.e.  any  two  points  in  M3  can  be  connected  by  a  single  geodesic  line.  Let  l(x,x 0) 
be  the  length  of  the  geodesic  line  connecting  points  x  and  xq.  Then  the  following  asymptotic 
behavior  of  the  function  w  and  its  derivatives  takes  place  for  |a|  <  2,  k  —  0, 1,  x  ^  xq 


DZDk,w(x,  s)  =  D°Dt 


exp  [—si  (. x ,  xo)] 

/  (x,  Xo) 


1  +  0  (  - 

s 


oo, 


(36) 


where  f  (. x ,  xq)  is  a  certain  function  and  f  (x,  xo)  ^  0  for  x  ^  xq. 

Having  the  data  at  only  one  side  Ti,  as  in  (8),  of  the  cube  is  not  sufficient  for  a 
good  stability  of  the  numerical  solution.  To  provide  a  better  stability,  we  now  derive  an 
approximate  boundary  condition  for  the  function  In  w  at  the  rest  T2  U  T3  of  the  boundary 
dii.  It  follows  from  (14)  and  (30)  that  the  function  w  satisfies  the  radiation  condition  at 
infinity,  lim^i^oo  (d\x\w  +  sw )  (x)  =  0,  where  d\x\ w  :=  drw  is  understood  in  terms  of  sperical 
coordinates  with  the  radius  r  :=  |x|  .  Hence,  assuming  that  the  number  P  in  (5)-(7)  is 
sufficiently  large  and  keeping  in  mind  that  we  work  with  an  approximate  mathematical 
model  (subsections  6.2,  6.3),  we  impose  the  following  approximate  boundary  condition  at 

r2  ur3 

( dnw  +  sw)  |r2ur3=  0.  (37) 

It  follows  from  (37)  that 

dn  (In  w  (x,  s))  |r2ur3=  ^s.  (38) 

Actually  condition  (38)  is  not  an  informative  one.  This  is  because  it  is  independent  on  the 
target  coefficient  c  (x) .  Hence,  it  is  logical  to  use  two  boundary  conditions  (8)  as  well  as  (38), 
since  this  should  likely  provide  both  a  better  reconstruction  accuracy  and  a  better  stability. 

Keeping  in  mind  our  target  application  of  imaging  of  antipersonnel  plastic  land  mines 
(section  1),  we  have  verified  the  approximate  boundary  condition  (37)  computationally,  both 
in  3-D  and  2-D  cases,  as  follows.  For  a  variety  of  cases  modeling  our  target  application 
(subsection  8.2)  we  have  computationally  solved  the  forward  problem  for  equation  (11)  in 
a  domain  H,  which  was  much  larger  than  the  domain  H  in  (5),  f!  C  Ll,dLl  fl  <91!  =  0. 
Because  of  (12),  we  have  imposed  the  zero  Dirichlct  boundary  condition  at  <912.  Next,  we 
have  solved  equation  (11)  in  the  domain  12  with  the  boundary  condition  (37)  at  T2  U  r3.  As 
to  Ti,  we  have  used  the  Dirichlet  boundary  condition,  which  was  calculated  from  the  above 


11 


solution  of  the  forward  problem  in  Q.  When  doing  so,  we  have  used  the  same  values  of  the 
parameter  s  for  which  we  have  numerically  solved  our  inverse  problem.  Comparison  of  these 
two  solutions  has  consistently  revealed  that  in  a  subdomain  Q  C  Q,  whose  boundary  had 
a  small  distance  from  T2  U  T3,  these  two  solutions  have  almost  coincided.  Thus,  the  above 
provides  a  numerical  justification  for  the  approximation  (37),  since  mine-like  inclusions  of 
our  interest  are  located  in  fi,  also  see  sub-subsection  8.1.1  for  a  relevant  discussion. 

4  The  Algorithm 

The  algorithm  of  this  section  was  described  in  detail  in  previous  publications  [4,  5]  for  the 
case  of  the  complete  data  collection  and  in  [30]  for  the  case  of  backscattering  data.  However, 
since  we  need  to  refer  to  some  parts  of  this  algorithm  in  our  convergence  analysis  in  sections 
5-7,  we  need  to  outline  it  briefly  here.  We  focus  only  on  those  items  which  we  need  in  sections 
5-7. 


4.1  The  sequence  of  elliptic  equations 

We  work  only  with  the  function  w(x,s),s  >  0  assuming  that  conditions  of  Theorem  3.1 
hold.  Since  by  (15)  w  >  0,  then  we  can  consider  the  function  v(x,s)  =  In  w/s2.  Let 
q  (. x ,  s)  =  dsv  (x,  s ) .  An  important  point  here  is  that 


v  (x,  s ) 


S 

fq(X,r)dr  +  V(X),V(X):=v(X,s). 

s 


(39) 


We  call  V  (x)  the  tail  function.  Here  s  >  1  is  the  truncation  pseudo  frequency  which  should 
be  chosen  in  numerical  experiments.  Actually,  s  is  one  of  regularization  parameters  of  our 
method.  Hence, 


V  (x)  :=  V  (x,  s)  =  s  2  In  w  (x,  s) .  (40) 

Assuming  that  conditions  of  Lemma  3.1,  we  obtain 

||H(a;,s)||c2+7(n)  =  O  (W1)  ,  IMIC2+7(n)  =  0  («~2)  , «  ^  °o-  (41) 

The  function  q  (x,  s )  satisfies  a  certain  nonlinear  integral  differential  PDE  in  which  Volterra- 
like  integrals  generated  by  (39)  are  involved.  We  have  one  equation  with  two  unknown 
functions  q  and  V.  The  reason  why  we  can  approximate  both  of  them  is  that  we  treat  them 
differently:  while  we  approximate  the  function  q  via  inner  iterations  “within”  that  equation, 
we  approximate  the  function  V  via  outer  iterations  solving  the  problem  (11),  (12)  and 
using  (40)  then.  Boundary  conditions  (8),  (38)  are  transformed  in  the  following  boundary 
conditions  for  q  (x,  s ) 

q  | rx  =  (x,  s ) ,  dzq  |ri=  i>i  (x,  s ) ,  dnq  |r2ur3=  s 


(42) 
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The  problem  of  approximating  functions  q  and  V  from  that  equation  is  the  major  dif¬ 
ficulty  of  our  method.  We  treat  it  via  a  layer  stripping  procedure  with  respect  to  s.  Let 
[s,  s]  C  (0,  cx))  be  the  s— interval  where  we  approximate  the  function  q  (. x ,  s ) .  In  practical 
computations  this  interval  should  be  chosen  via  numerical  experiments  conducted  for  a  range 
of  parameters,  which  is  suitable  for  a  specific  application  of  ones  interest.  We  approximate 
the  function  q  (x,  s )  as  a  piecewise  constant  function  with  respect  to  the  pseudo  frequency 
s  G  [s,  s].  That  is,  we  assume  that  there  exists  a  partition  s  =  sn<  sjv-i  <  ...  <  si  <  sq  =  s 
of  the  interval  [s,  s]  with  the  sufficiently  small  grid  step  size  h  =  s*_i  —  st  such  that 
q(x,s )  =  qn(x)  for  .s  G  (sn,sn„i],n  =  1, ...,  N.  Next,  we  multiply  the  above  mentioned 
equation  for  q  by  the  s-dependent  CWF  Cn tfl  (s)  =  exp  [—fj,  |s  —  sn_i|]  ,  s  G  (sn,  sn_i]  ,  /i  >>  1 
and  integrate  with  respect  to  s  G  (sn,  sn_i] .  The  parameter  //  should  be  chosen  in  numer¬ 
ical  experiments.  As  a  result,  we  obtain  a  coupled  system  of  elliptic  PDEs  for  functions 
qn.  Because  of  Volterra-like  s— integrals  in  the  original  equation  for  q,  this  system  can  be 
solved  sequentially  starting  from  q1.  When  solving  these  equations,  we  have  m  >  1  inner 
iterations  for  each  n  with  respect  to  the  tail  functions.  Thus,  for  each  n  G  [1,  iV]  we  obtain 
a  finite  sequence  of  functions  qUtk,  k  G  [1,  m] .  Those  equations  and  boundary  conditions 
for  functions  qn^  are 


A  q. 


n,k 


I  V Qn,k  —  ^2 


71—1 


+2i42>„VVrn,jfe  ^  “  A2 ,n  (yvn,k)2  ,X  G  (n,  k)  G  [1,  N]  x  [1,  m\ , 

9n,fc|ri  V^0,7i(*^)  7  &zQn,  A;  |  F  i  *01,71  (*^)j  9nQn,k  |r2UP3  (^n^n— l) 


(43) 

(44) 


Functions  and  are  averages  over  the  interval  (sn,sn_i)  of  respectively  func¬ 

tions  ip o  (x,  s )  and  ipi  (x,  s )  in  (42)  and  are  therefore  approximate  boundary  conditions.  The 
third  boundary  condition  in  (44)  is  the  average  of  the  function  s~2  in  (42)  over  the  interval 
(sn,  sn- 1) .  In  (43)  Ai;n  =  Ai  n  (. h ,  /x)  >  0,  i  =  1,  2  are  known  parameters,  see  specific  formulas 
in  [4].  It  is  convenient  to  set  in  (43) 

q0  =  0.  (45) 


Remarks  4.1: 

1.  In  fact,  one  should  also  have  the  nonlinear  term  Bn  (/i,  h )  (' Vqn,k )2  in  the  right  hand  side 
of  equation  (43)  with  a  certain  parameter  Bn  ( n,h )  [4],  However,  because  of  the  presence 
of  the  CWF  we  have  \Bn(fj,,h)\  <  8s2//i  for  fill  >  1  [4],  We  have  used  fi  =  50  in  our 
computations.  Hence,  assuming  that  yU  >>  1,  we  ignore  this  term.  This  allows  us  to  solve 
a  linear  problem  for  each  qnj,.  We  have  also  conducted  numerical  experiments  for  the  case 
when  this  term  is  not  ignored.  In  this  case  we  have  also  performed  iterations  with  respect 
to  this  term.  However,  results  of  those  numerical  studies  have  shown  that  the  influence  of 
this  term  is  negligible  for  yU  >  50. 
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2.  In  principle  one  can  avoid  using  the  CWF  via  choosing  the  size  h  of  each  subin¬ 
terval  (sn,  sn-i)  to  be  very  small.  This  would  again  result  in  a  small  multiplier  at  (V^)2. 
However  the  number  N  of  subintervals  would  increase  in  this  case,  which  would  make  it 
necessary  to  solve  more  problems  for  functions  qn.  Hence,  both  the  computational  time  and 
the  computational  error  would  increase. 

3.  Setting  Bn  i^7qn,k)2  :=  0  for  large  values  of  //  is  not  a  linearization.  The  nonlinear¬ 
ity  actually  surfaces  in  the  iterative  process  in  those  terms  of  (43)  which  contain  q3  and 
Vn>k-  Indeed,  these  terms  contain  products  Vq3Vqn,k-  Also,  tails  Vn^  depend  nonlinearly  on 
functions  qj,qn,k ,  see  (50)  and  (51). 

4.  In  our  analytical  derivation  we  use  the  third  condition  (42)  at  T3  only  at  s  —  s,  where 
s  is  the  highest  value  of  s  we  use  in  our  algorithm,  see  (104),  (105)  and  Lemma  7.1.  However, 
we  use  the  third  condition  (42)  at  T2  for  all  values  of  s  in  our  convergence  analysis.  As 
to  the  numerical  studies,  our  computational  experience  has  shown  that  it  is  better  for  the 
stability  to  use  the  third  condition  (42)  at  T2  U  T3  for  all  values  of  s,  and  this  is  what  we  do 
in  our  numerical  testing,  see  subsection  8.1  for  a  relevant  discussion. 

5.  In  addition,  if  using  the  Laplace  transform  of  only  one  of  conditions  (8)  combined 
with  the  third  condition  (42),  then  there  is  an  analytical  difficulty  here,  compared  with  [4,  5]. 
Indeed,  in  the  convergence  analysis  of  [4,  5]  the  boundary  <90  G  C3,  which  enabled  one  to  use 
Schauder  theorem.  However,  the  boundary  of  our  cube  O  in  (5)-(7)  is  not  smooth.  Hence, 
the  existence  of  a  smooth  solution  for  the  corresponding  elliptic  boundary  value  problem  for 
the  function  in  section  4  cannot  be  guaranteed  in  this  case. 

6.  In  summary,  it  follows  from  the  above  two  remarks  that  it  is  optimal  in  computations 
to  use  the  Laplace  transforms  of  both  boundary  conditions  (8)  as  well  as  the  third  condition 
(42).  The  over- determination  in  boundary  conditions  (8)  logically  leads  to  the  QRM. 

The  iterative  process  of  the  next  subsection  reconstructs  iterative  approximations  cn^.  (x)  G 
C 7  (O)  of  the  function  c  (x)  and  cnjfc  (x)  G  [1,  d]  in  O.  On  the  other  hand,  to  iterate  with  re¬ 
spect  to  the  tails,  we  need  to  solve  the  forward  problem  (11),  (12)  on  each  iterative  step.  To  do 
this,  we  extend  each  function  cn^  (x)  outside  of  the  domain  O.  So  that  the  resulting  function 
cn,fc  (x)  =  1  outside  of  O,  cn)fc  (x)  =  cn ^  (x)  in  a  subdomain  O'  C  O,  dntk  (x)  G  [1,  d]  ,  Wx  G  M3 
and  crhk  G  C7(M3).  When  estimating  solutions  resulting  from  the  QRM,  we  need  to  use 
a  Carleman  estimate,  which,  however,  estimates  the  solution  only  in  a  subdomain  of  the 
original  domain.  The  way  out  of  this  is  that  in  applications  one  can  always  choose  such  a 
domain  of  interest  which  is  a  little  bit  bigger  than  the  one  where  the  coefficient  is  unknown, 
and  we  need  this  in  our  convergence  analysis  in  section  7.  These  considerations  prompt  us 
to  apply  the  construction  described  below  in  this  subsection. 

Let  P  be  the  number  in  (5)-(7)  and  P0  =  const.  G  (0,  P ) .  Denote  Op0  =  On{z  G  (0,  P0)}  • 
We  will  assume  below  that  c  (x)  =  1,  Wx  G  M3\Op0.  Consider  a  subdomain  O'  C  Op0.  Choose 
a  function  xi  (x)  G  C°°  (1R3)  such  that 


f  1  in  O', 

Xi  ( x )  =  <  between  0  and  1  in  Op0\O', 
[  0  outside  of  Op0. 


14 


The  existence  of  such  functions  \i  (x)  is  well  known  from  the  Real  Analysis  course.  Define 
the  target  extension  of  the  function  critk  as 

Cn,k(x)  '■=  (1  —  Xi  (x))  +  Xi  (x)  cn,k  (x) ,  Vx  G  M3.  Hence,  cUik(x)  =  1  outside  of  QPo  and 
cn,fc  G  C1  (M3).  Using  cUtk  ( x )  G  [1,  d] ,  Wx  G  H,  it  is  easy  to  verify  that  cn tk  (x)  G  [1,  d]  ,  Vx  G 
M3. 

4.2  The  iterative  process 

First,  we  choose  an  initial  tail  function  V\,\  ( x )  G  C2+7  (fl).  Onr  choice  of  this  function  is 
described  in  subsection  6.3  and  it  is  different  from  ones  used  in  previous  publications  on 
this  method.  Let  m  >  1  be  an  integer  which  we  choose  in  numerical  experiments.  For 
each  n  G  [1,  IV]  we  have  m  inner  iterations  with  respect  to  the  tails  via  computing  functions 
qn,k,Vn,k,k  G  [1,  m]. 

Step  nk ,  where  n  G  [1,  N]  ,  k  G  [1,  m]  .  Recall  that  by  (45)  q0  =  0.  Suppose  that  functions 
q.j  G  H5  (fl)  ,j  =  l,...,n  —  1  and  tails  Vi, Vn-i,  Vn^  G  C2+7  (fl)  are  constructed.  First, 
we  change  the  operator  and  the  right  hand  side  of  (43)  as  described  in  subsection  4.3.  To 
construct  the  function  qn.k,  we  use  the  QRM  (subsection  4.3)  to  find  the  “least  squares” 
solution  of  the  perturbed  boundary  value  problem  (43),  (44).  Hence,  we  obtain  the  function 
qn^k  G  H5  (f2) .  To  reconstruct  an  approximation  cU)k  ( x )  for  the  function  c  ( x ) ,  we  first  use 
our  discrete  analog  of  (39)  to  calculate  an  approximation  for  the  function  v  ( x ,  sn)  and  then 
calculate  an  approximation  cn^  for  the  function  c.  Specifically, 

n— 1 

Vn,k  {x,  Sn )  =  -hqHik  (x)  -hj^  qj  (x)  +  Vn>k  (x) ,  x  G  flPo,  (46) 

3= 0 

Cn,k  (*^)  An n  k  (^r,  Sn)  T  sn  |  V vn^k  (x,  sn)  |  ,  x  G  .  (47) 

Since  we  assume  in  (3)  that  the  exact  solution  of  our  Inverse  Problem  c*  G  [1,  d]  and  since 
the  computed  function  cn,fc  in  (47)  does  not  necessarily  satisfy  this  requirement,  we  set 

{Cn,k  (x)  ,  if  Cn,fc  (x)  G  [1,  d] ,  x  G  Dp0, 

1,  if  cn,fc  (x)  <  1,  x  G  DPo,  (48) 

d,  if  cn,fc  (x)  >  d,x  G  flp0. 

Since  functions  qj,  qUtk  G  i/5  (D)  ,  then  the  embedding  theorem  implies  that  qv  qnP  G  U3  (fl)  . 
In  addition,  the  tail  function  VU)k  G  C2+7  (O)  .  Hence,  (46)- (48)  imply  that 

Cn,k  e  C7  (H)  .  (49) 

Next,  we  construct  the  function  dn,tk  (x)  as  explained  in  subsection  4.1.  By  (49)  the  function 
Cn^eF7  (M3) .  Next,  we  solve  the  forward  problem  (11),  (12)  with  c  (x)  :=  cntk  (x)  for  s  :=s 
and  obtain  the  function  wnjk  (x,  s)  .  Existence  and  uniqueness  of  the  function  wn^k  (x,  s )  in 
the  form  (13),  (14)  are  guaranteed  by  Theorem  3.1.  We  set  for  the  new  tail 

VHtk+ 1  (x)  =  s~2  In  wntk  (x,  s)  G  C'2+1  (fl)  if  k  <  m. 


(50) 
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And  if  k  —  m ,  then  we  set  cn  ( x )  :=  cn.m  (x)  ,  qn  (. x )  :=  qn>m  (x).  And  for  the  tail 

Vn  {x)  :=  K,m+i  (x)  :=  s~2  In  wn,ra  (x,  s)  :=  K+1,1  (x)  for  x  e  Cl.  (51) 


4.3  The  quasi-reversibility  method 

Since  the  Carleman  estimate  of  section  5  estimates  the  target  function  only  in  a  subdomain 
of  the  domain  fi,  we  introduce  the  function  \2  (x)  which  is  the  characteristic  function  of  the 
subdomain  Clp0,X2  (x)  =  1  in  Qp0  and  \2  (x)  =  0  in  IR3\flp0.  Denote 


®n,fc  (x) 


Hn,k  (x) 


n—  1 


M,n  X2  (x)  h  ^  -  Vbi.lt  , 


3=0 


^  n—  1 


n—1 


-v2  Ev%,  X2  (x)  +  2 A2,nVVn  k  /j  V(lj  X2 


X) 


0=0 


3=0 


~A2,n  (Wn ,kf 


(52) 

(53) 


Hence,  the  function  Hn^{x)  is  the  perturbed  right  hand  side  of  equation  (43).  Note  that 
Hn.k  E  L2  (Cl) .  The  perturbed  boundary  value  problem  (43),  (44)  is 


^ Qn,k  ®"n,k^  Qn,k  Hn,ki  (54) 

Qn.k  | F]  V’O ,n(x)j  9zqn,k  | F i  f/^1  ,n(x)j  9nqn  k\v‘2\jTz  (^n^n—l)  •  (55) 

Since  we  have  two  boundary  conditions  rather  then  one  at  lb,  we  find  the  “least  squares” 
solution  of  the  problem  (54),  (55)  via  the  QRM.  Specifically,  we  minimize  the  following 
Tikhonov  functional 


Jn,k(u)  =  II- Au  -  an,kX7u  -  Hn,k\\lm  +  a  \\u\\2H5m  ,  (56) 

subject  to  the  boundary  conditions  (55).  Here  a  G  (0, 1)  is  the  small  regularization  parame¬ 
ter.  Let  u  (x)  be  the  unique  minimizer  of  this  functional,  which  is  guaranteed  by  Lemma  5.2. 
Then  we  set  qnj .  (x)  :=  u  (x)  .  Local  minima  do  not  occur  here  since  (56)  is  the  sum  of  square 
norms  of  two  expressions,  both  of  which  are  linear  with  respect  to  u.  The  second  term  in 
the  right  hand  side  of  (56)  is  the  Tikhonov  regularization  term.  We  use  the  H 5  (0)  —norm 
here  in  order  to  ensure  that  the  minimizer  u  :=  qn y-  G  C 3  (0) .  It  was  shown  above  that  the 
latter  implies  (49).  We  call  the  minimizer  u(x)  of  the  functional  J“fc(u)  the  QRM  solution 
of  the  problem  (54),  (55). 

5  Estimates  for  the  QRM 

For  brevity  we  scale  variables  in  such  a  way  that  in  (5)-(7)  P  =  1/2  in  sections  5-7 

Cl  =  {x  =  (x,  y,  z)  :  (x,  y )  G  (-1/4, 1/4)  x  (-1/4, 1/4) ,  z  G  (0, 1/2)}  (57) 
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Let  A,  v  >  2  be  two  parameters.  Introduce  the  z— dependent  CWF  K(z), 

K  (z)  :=  K\jV(z)  =  exp(A p_I/),  where  p(z)  =  z  +  1/4,  z  >  0. 


This  CWF  is  different  from  the  ones  previously  used  for  Carleman  estimates  for  elliptic  PDEs 
[16,  28,  34],  Note  that  p  (z)  G  (0,  3/4)  in  Q  and  p  (z)  |p3=  3/4.  Let  the  number  x  G  (1/3, 1) . 
Denote  =  {x  E  Q  :  p(z)  <  3x/4} .  Hence,  if  xi  <  x2,  then  C  Also,  Hi  =  D  and 
fli/3  =  0.  In  addition,  K2  (z)  >  exp  [2A  ((4/3)  x) V(: }  in  Everywhere  below  C  >  0  denotes 
different  generic  constants  which  are  independent  on  the  domain  0  in  (57),  (•,  •)  denotes  the 
scalar  product  in  L2  (O)  and  [•,  •]  denotes  the  scalar  product  in  H 5  (D). 

Lemma  5.1.  Fix  a  number  v  :=  u0  >  2.  Consider  functions  u  G  H3  (D)  such  that  ( see 

m  v) 

U  |ri=  Uz  |ry  =  dnu  |r2=  0.  (58) 

Then  there  exists  a  constant  C  >  0  such  that  for  all  A  >  2  the  following  Carleman  estimate 
is  valid  for  all  these  functions 

f  (Au)2K2dx  >  j  [  (D“m)2  1<2(ix  +  c  [  [A  (V«)2  +  A3m2]  K2dx 
n  A  M=2  n  n 

-C\3  ||u||H3(n)  exP  [2A  (4/3 H  • 

Proof.  It  is  convenient  to  initially  assume  that  v  >  2  is  an  arbitrary  number.  We  have 


(Am)2  K 2  =  (u2xx  A  u2yy  A  u2zz  A  2uxxuzz  +  2 uxxuyy  A  2uyyuzz)  K2 

=  {UL  +  uly  +  ulz)  I{2 

+dx  (2 uxuzzK2  +  2 uxuyyK2)  A  <%,  (2 uyuzzK2) 

2nxUZZXI\.  2 UXUyyXI\.  ^  2UyVjZZyI\. 

=  (u2xx  A  u2y  A  m22)  A'2  +  dx  (2uxuzzK2  +  2uxuyyK2)  +  dy  (2uyuzzK2) 
+dy  (~2uxuxyK2)  +  2u2xyK2 
+dz  (- 2uxuxzK 2)  +  2uxzK2  -  4:\up~u~luxuxzK 2 
AcL  (— 2uyuyzK 2)  +  2uyzK2  —  A\vp~u~luyuyzK2 . 


Thus,  we  have  obtained  that 

(Am)  A 2  =  (m2^  +  M2y  +  m22  +  2u2y  +  2m2„  +  2m2.,)  K2 

4Az/p  ( 'U'x'Ujxz  A  tLyiiyz )  A  T  c/e  [2  {uxuzz  T  uxuyy )  A  J  (59) 

A ciy  [2  ( uyuzz  ‘Ujx'u,xy)  A +  cb  [  2  (m xuxz  A  UyUyzf  K  j  . 

Using  the  Cauchy-Schwarz  inequality  “VA  >  0”,2a6  >  — ea2  —  b2  /e  and  taking  £  =  1/4,  we 
obtain 


4Az/p  "  1  {uxuxz  +  uyuyz)  K 2  >  -  (m22  +  m2J  A'2  -  4AVp  2l/  2  (Vm)2  A'2. 
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Hence,  (59)  implies  that 

(Am)2  K2  >  ^  ( Dauf  K2  -  4A  Vp"2^2  (Vm)2  K2  +  dx  [. 2ux  (uzz  +  mw)  A'2] 


|a|=2 

~\~dy  [2  ( UyUZZ  Ux^xy)  A  ^  2  ( 'U'x'U'XZ  A  MyMy_2 )  A"  j  .  (60) 

Consider  a  new  function  v  =  mAT.  Substituting  u  =  vK~l ,  we  obtain 

(Am) V+1A2  =  (yi  +  y2  +  y3)V+1  >  2y2(yi  +  p;^1,  (61) 

Pi  =  Am,  p2  =  2Xvp~l/~1vz,  y3  =  (Az/)2p“2l/_2(l  -  (1/  +  1)  (Az/)_1  /A)u.  (62) 

We  have 

2pi y2pu+1  =  dx  (AXvvzvx)  +  dy  (■ AXvvzvy )  +  &  [2Ai/  (u2  -  u2  -  u,2)]  .  (63) 

Next,  by  (61)  and  (62) 

2y2y3pv+1  =  A(\v)3  (p-2"-2  -  (v  +  1)  (Xv)-1  p-"-2)  vzv 

=  dz  [2(Az/)3  ( p~2v~ 2  ~{y+  1)  (Az/)_1  p-*'-2)  u2]  (64) 

+4(Az/)3  (1/  +  1)  p~2u~ 3  (1  -  (1/  +  2)  (2Az/)_1  p")  u2 
>  2A3z/4p_2l/_3M2  +  [2(Az/)3  (p"2y“2  -  (1/  +  1)  (Az/)_1  p"^2)  m2]  . 

Summing  up  (63)  with  (64),  using  (61)  and  the  backwards  substitution  u  =  vK,  we  obtain 

(A ufpv+1K2  >  2X3u4p-2v~3u2K2  +  dxUx  +  3yA2  +  dzU3,  (65) 

where  the  following  estimates  are  valid  for  functions  U \ ,  f/2,  U3 

\Ui\  <  CXv  \ux\  (\uz\  +  Xvp~v^  |m|)  A'2, 

|A2|  <  CXv  \uy\  (|mz|  +  Xvp~u~1  |m|)  A'2,  (66) 

\U3\  <  CXv  (|Vm|2  +  X2u2p~2v-2u2)  K2. 


We  now  need  to  incorporate  the  term  A  (Vm)2  K 2  in  the  right  hand  side  of  the  Carleman 
estimate.  Hence,  we  continue  as  follows: 

— XvuAuK2  =  dx  (— XvuuxK2)  +  dy  (—XvuuyK2)  +  dz  (— XvuuzK2) 

+Xv  (Vm)2  K 2  -  2X2v2 p~u~1uzuK2  =  Xv  (Vm)2  K2 

-2A3z /3p-2"-2  (1  +  {v  +  1)  (2Az/)-1  pv)  u2K 2  +  dxUA  +  ayV5  +  dzU6, 

Hence, 

- XvuAuK 2  >  Xv  (Vm)2  A'2  -  4A3z/3p-2"-2M2A'2  +  dxU4  +  dyU5  +  dzU6,  (67) 
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IJ4  =  -X vuuxK2,  U5  =  - \vuuyK 2,  \U6\  <  C  {Xvul  +  AV/T^V)  K 2.  (68) 

Sum  up  (65)  and  (67),  take  into  account  (66)  and  (68)  as  well  as  the  fact  that 

2A 3v4p~2v~3  -  4A 3v3p~2u~2  =  2A 3v4p~2u~3  (l  -  p  (2 z/)_1)  >  A 3v4p~2v~3. 

We  obtain 

(A  u)2K2  -  XvuAuK 2  >  (69) 

Xu  (Vw)2  K 2  +  X3  v4  p~2u~3u2  K2  +  dxU7  +  dyU8  +  dzU9, 

\U7\  <  CXv  \ux\  ( \uz\  +  A vp~u~l  |u|)  K2,  (70) 

\U8\  <  CAz/  l-u.yl  (|uz|  +  A vp~v~l  |m|)  K2 \  (71) 

\Ug\  <  CXv  (|  Vu|2  +  A 2v2p-2u~2u2)  K2.  (72) 

Fix  the  number  v  :=  u9  >  2.  Then  we  can  incorporate  u9  in  C  and  can  regard  that  puo+1  < 
C ,  since  puo+1  <  1.  By  the  Cauchy-Schwarz  inequality  —XvuAuK2  <  (Au)2 pu+1  K 2/2  + 
X2v2 p~u~1u2K 2/2.  Hence,  we  obtain  from  (69)-(72) 

(A  ufK2  >C[  X  (V«)2  +  A  V]  K2  +  4H7  +  dyU8  +  dzUg.  (73) 

We  now  set  in  (60)  v  v0,  divide  it  by  A p  with  a  positive  constant  p  =  p[y 0)  such  that 
4vq p-2v°~2 /m  <  (7/2,  add  the  resulting  inequality  to  (73)  and  take  into  account  (70)-(72). 
Then  with  a  new  constant  C  we  obtain  the  following  pointwise  Carleman  estimate  for  the 
Laplace  operator  in  the  domain  72 

(A  ufK2  >  j  J2  ( Dauf  K 2  +  C  [X  (V«)2  +  A3u2]  K 2  (74) 

\a\=2 

+dxU\  0  +  dyUn  +  dzUi2, 

| TT 10 1  <  CX  | ux\  (|uzz|  +  \uyy\  +  \uz\  +  A  \u\)  K2 ,  (75) 

|Hn|  <  C  [A  \uy\  (\uzz\  +  \uz\  +  A  \u\)  +  \uxy\  \ux\]  K2  (76) 

\U12\  <  CX[\uxz\2  +  \uyz\2+\Vu\2  +  X2u2}K2.  (77) 

We  now  integrate  the  formula  (74)  over  the  rectangle  72  using  the  Gauss  formula.  It  is 
important  that  because  of  (58)  and  estimates  (75)-(77),  resulting  boundary  integrals  over  Ti 
and  T 2  will  be  zero.  We  obtain 

I {Au)2K2dx  >  j  f  (Dauf  K2dx  +  C  J  [X  (V«)2  +  A3u2]  K2dx  (78) 
n  l“l=2  n  n 

-CX  I  [\uxz\2  +  \uyz\2  +  |Vn|2  +  AV]  K2dS. 
r3 

Note  that  K2  (1/2)  =  K 2  (. z )  |r3=  riling  K2  (. z )  =  exp  [2A  (4/3)^°]  .  Thus, 

J  X  +  \uyZ | 2  +  | X7u\2  +  A2u2]  K 2dx  <  CX3  ||u||^-3^n^  exp  [2A  (4/3 )i/°] . 
r3 
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Substituting  this  in  (78),  we  obtain  the  estimate  of  this  lemma.  □ 

A  peculiarity  of  the  proof  of  Lemma  5.1  is  that  when  integrating  the  pointwise  Carleman 
estimate  (74)  over  the  domain  fl,  we  should  take  into  account  only  one  (rather  than  conven¬ 
tional  two)  zero  boundary  condition  (58)  at  T2.  This  requires  a  careful  evaluation  of  terms 
under  signs  of  derivatives  dx,  dy  which  was  not  done  before. 

We  now  establish  both  existence  and  uniqueness  of  the  minimizer  of  the  functional  (56). 
Denote  a^\  (x)  ,i  —  1,2,3  components  of  the  vector  function  aUjk  (x)  in  (52),  (54).  Let 

a:)),  <  M,  M  =  const.  >  0,  %  —  1, 2, 3.  (79) 

’  noo(n) 

Lemma  5.2.  Assume  that  there  exists  a  function  G  H5  (f2)  satisfying  boundary  condi¬ 
tions  (55),  except  of  maybe  at  T3,  and  condition  (79)  holds.  Then  there  exists  unique  mini¬ 
mizer  u  G  H5  (D)  of  the  functional  (56).  Furthermore,  with  a  constant  C±  =  C\  (M)  >  0 

IMItf5(n)  —  C\Oi  ^  (|I^4'IIl2(0)  +  ll^>llrr5(n))  • 

Proof.  Let  u  be  a  minimizer  of  J£k(u)  satisfying  boundary  conditions  (55).  Let  U  = 
u  —  <L.  Then  the  function  U  satisfies  boundary  conditions  (58).  By  the  variational  principle 

( Gn,kU ,  Gn)kv)  +  a  [U,  v ]  =  (Hnik  -  Gn,k§,  Gn,kv)  -  a  [<L,  v] , 
for  all  functions  v  G  H 5  (D)  satisfying  boundary  conditions  (58).  Here 

Gn,kU  :=  AU  —  an^kSjU.  (80) 


The  rest  of  the  proof  follows  from  the  Riesz  theorem.  □ 

Lemma  5.3  ([30],  see  the  second  paragraph  of  Introduction) .  Let  the  function  u  G  H 5  (O) 
satisfies  boundary  conditions  (58)  as  well  as  the  variational  equality 

0 Gn)ku ,  G,hkv)  +  a  [u,  v ]  =  {HnM,  Gn,kv)  +  a  \g,  v]  (81) 


for  all  functions  v  G  H5  (D)  satisfying  (58),  where  the  operator  Gn  k 
Then 


\u\\h5(q.)  —  a  1 ^ 


*n,k  II l2(0)  +  \\9\\h5(Q.)  ■ 


is  defined  in  (80). 


Theorem  5.1  ([30],  see  the  second  paragraph  of  Introduction).  Assume  that  condition 
(79)  holds.  Let  g  G  H5  (D)  be  an  arbitrary  function.  Let  u  G  Hb  (O)  be  the  function 
satisfying  boundary  conditions  (58)  as  well  as  the  variational  equality  (81)  for  all  functions 
v  G  H5  (f2)  satisfying  (58).  Let  the  number  x  G  (1/3, 1)  and  the  number  p  G  (x,  1) .  Consider 
the  numbers  bi,b2, 


(82) 
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where  u0  is  the  parameter  of  Lemma  5.1.  Then  there  exists  a  sufficiently  small  number 
a i  =  ai(M,  x,  p)  G  (0,1)  such  that  for  all  a  G  (0,  af)  the  following  estimate  holds  with  a 
constant  C2  =  C2  (M,  x,  p)  >  0 

IMI/rqn*)  —  1  || l2(h)  +  a  2  llfi,lli/5(n)  •  (88) 

The  proof  of  Theorem  5.1  is  based  on  Lemma  5.1.  An  important  difference  between 
Theorem  2.1  of  [16]  and  Theorem  5.1  is  that  we  now  use  a  new  CWF  in  order  to  take 
into  account  the  Neumann  boundary  condition  at  T2-  Indeed,  an  analog  of  our  Neumann 
boundary  condition  at  T2  was  not  used  in  [16].  Recall  that  we  need  the  Neumann  boundary 
condition  at  T2  only  to  ensure  a  better  stability  of  the  numerical  solution  of  our  target  CIP 
(the  fourth  Remark  4.1). 

6  Preliminaries  for  the  Convergence  Analysis 

We  follow  the  concept  of  Tikhonov  for  ill-posed  problems  [19,  34,  40].  By  this  concept,  one 
should  assume  that  there  exists  an  “ideal”  exact  solution  of  an  ill-posed  problem  with  the 
“ideal”  exact  data.  Next,  one  should  prove  that  the  regularized  solution  is  close  to  the  exact 
one.  In  the  course  of  this  proof  one  can  use  an  a  priori  upper  estimate  of  a  certain  norm  of 
the  exact  solution. 

6.1  Exact  solution 

First,  we  need  to  introduce  the  definition  of  the  exact  solution.  Several  aspects  of  this 
definition  are  different  from  ones  in  previous  publications  [4,  5,  27].  We  assume  that  there 
exists  a  coefficient  c*  ( x )  which  is  the  unique  exact  solution  of  our  Inverse  Problem  with  the 
exact  (i.e.  noiseless)  data  Tp*0  (x,t) ,  </?*  (x,t)  in  (8)  (Remark  2.1).  We  assume  without  loss  of 
generality  that 

llc*llcu(R3)  —  ^  and  conditions  (3),  (4)  hold  for  c*  (x) .  (84) 

Let  u*  ( x ,  t)  be  the  solution  of  the  forward  problem  (1),  (2)  with  c  :=  c*  and  w*  (. x ,  s)  be  its 
Laplace  transform  (10)  for  s  >  s°  (d)  >  0.  Since  the  source  Xq  fi  Q,  then  it  follows  from 
(84)  and  Theorem  3.1  that  w*(x,s)  G  C5+1  (hi)  .  Similarly  with  section  4  let  v*(x,s)  = 
s~2  In  w*  (x,  s )  and  q*  (. x ,  s)  =  dsv*  (x,  s) ,  s  >  s(d) .  Let  [s,  s]  be  the  s— interval  of  section  4 
and  s>s(d).  We  assume  that 

q*  g  C5+7  (n)  X  C1  [s,s] ,  ||g*  (x,s)||c5+7(n)xci[1,3]  ^  c**  c*  =  const •  >  °>  (85) 

where  C*  is  the  given  upper  bound  of  the  norm  in  (85).  Consider  the  same  partition  of 
the  interval  [s,  s]  into  N  small  subintervals  as  in  section  4.  Let  q*  (. x )  be  the  average  of  the 
function  q*  (x,  s )  over  the  interval  (sn,  sn- 1)  .  Then  (85)  implies  that 


max  ||  q*(x)~  q*  (x,  s)||c5+7m)  <  C*h ■ 

SG[sn,Sn-lJ  V  / 


(86) 
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Hence,  increasing,  if  necessary  the  number  C*,  we  can  assume  that 

<  C".  (87) 

Let 

$5  (u«)  =  <1*  (x>s)  ki,  V’*  (x,s)  =  dzq*  (x,s)  |ri,  se[s,s].  (88) 

For  x  G  Tx  let  functions  (x)  and  ip*  n  (. x )  be  averages  of  functions  ipQ  {xi  s )  and  t/’i  0*5  s) 
over  the  interval  (sn,  sn_i)  .  Then  boundary  conditions  for  functions  g*  (x)  at  Ti  are 

?nlr!  =  ^0 ^Snlri  =  (89) 

The  exact  tail  function  V*  ( x )  is 

V*  (x)  =  s~2  In  w*  (x,  s) .  (90) 

The  function  q*  satisfies  the  following  analogue  of  equation  (42) 


71—1 


'  n—  1 


A q*n  -  Alin  h  Vg;  (x)  -  VU*  Vq*n  =  -A2,nh2  Vg* 


x) 


3=0 
n—  1 


W=0 


+  2  A2,nVV*  h  Ev«;  (x)  )  -  A2,n  (VU*)2  +  Fi,n  (x,  h,  fi) ,  q*  :=  0. 

V  1=0 

Similarly  with  (46)  and  (47) 


71—1 


<  M  :  =  ~hq‘n  (x)  -  h  ^  q]  (x)  +  V’  (z) ,  x  €  Si, 

c*  (x)  =  Ax*  (z)  +  s);  Vx*  (x)|2  +  F2,n  (x,  h)  ,x  e  U. 


(91) 


(92) 

(93) 


In  (91)  and  (93)  functions  F\ )Tl  (x,  h,  /i) ,  F2;n  (x,  h)  represent  approximation  errors  due  to  av¬ 
eraging  and  the  s— dependent  Carleman  Weight  Function.  In  particular,  the  term  Bn  (/i,  h)  (Vg*)2  , 
which  we  have  ignored  in  (43)  (Remarks  4.1),  is  a  part  of  Fj>n.  Although  we  can  prove  an 
analog  of  Theorem  7.1  for  the  case  F\  n  ^  0,  F2n  ^  0,^Qn  7^  ip0tn,ipln  7^  t/h,n,  this  would 
require  more  space.  At  the  same  time,  the  method  of  the  proof  would  be  almost  the  same. 

Thus,  we  will  “allow”  the  error  in  the  boundary  data  at  Ti  to  be  present  only  at  s  :=  s,  see 
Lemma  7.1.  Therefore,  for  brevity  only  we  assume  below  that 


F\,n  =  0,  F2,n  =  0,  -00, n  =  1p0 ,n,  V'gn  =  ^1  ,n,  ™  £  [1,  A]  .  (94) 

Using,  ideas  of  the  proof  of  Theorem  3.1,  it  is  possible  to  prove  that  for  s  >  s(d)  not 
only  the  function  w*  (x,  s )  G  U5+7  (fl)  but  also  functions  Dksw*  (x,  s )  G  U5+7  (fl)  ,  k  —  1,  2. 
Since  this  implies  that  g*  (x,  s)  G  C5+7  (fl)  x  Cl  [s,  s\ ,  then  it  is  not  necessary  to  impose  this 
assumption  in  (85).  However,  we  still  prefer  to  use  this  assumption  in  order  to  save  space 
and  also  because  this  is  not  our  main  focus.  The  reason  why  we  require  the  C4—  smoothness 
of  c*  in  (4)  is  to  ensure  that  V*  G  C5+1  (fl)  .  We  need  the  latter  to  justify  the  assumption 
in  (97)  that  p*  G  H5  (Q). 
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6.2  Approximate  global  convergence 

Definition  6.1  (approximate  global  convergence).  Consider  a  nonlinear  ill-posed  problem. 
Suppose  that  a  certain  approximate  mathematical  model  M  is  proposed  to  solve  this  problem 
numerically.  Assume  that,  within  the  framework  of  the  model  M,  this  problem  has  unique 
exact  solution  x*  €  B  for  the  noiseless  data  y*.  Here  B  is  an  appropriate  Banach  space  with 
the  norm  ||-||B .  Consider  an  iterative  numerical  method  for  solving  that  problem.  Suppose 
that  this  method  produces  a  sequence  of  points  {xn}^=1  C  B ,  K  G  [1,  oo) .  Let  the  number 
6  €  (0, 1) .  We  call  this  numerical  method  approximately  globally  convergent  of  the  level  6  in 
the  norm  of  the  space  B,  or  shortly  globally  convergent,  if,  within  the  framework  of  the  model 
M,  a  theorem  is  proven,  which  claims  that,  without  any  knowledge  of  a  sufficiently  small 
neighborhood  of  x*,  there  exists  a  number  K  e  [1,  K)  such  that  the  following  inequality  is 
valid 

11b  <6,Vn>K.  (95) 

1 1  ^ n  1 1 B 

Suppose  that  iterations  are  stopped  at  a  certain  number  k  >  K.  Then  the  point  Xk  is  denoted 
as  Xk  Xgiob  and  is  called  the  approximate  solution  resulting  from  this  method. 

This  is  our  formal  mathematical  definition  of  the  approximate  global  convergence  prop¬ 
erty.  However,  since  the  approximate  mathematical  model  M  is  involved  in  it,  then  it  is 
reasonable  to  raise  the  question  about  a  verification  of  M.  It  is  natural  to  address  this 
question  via  computational  testing.  Hence,  the  following  two  informal  conditions  should  be 
added  to  Definition  6.1: 

A.  Numerical  studies  should  confirm  that  xgiQb  is  indeed  a  sufficiently  good  approximation 
for  the  exact  solution  x* . 

B  (optional).  Testing  of  this  numerical  method  on  appropriate  experimental  data  also 
demonstrates  that  iterative  solutions  provide  a  good  approximation  for  the  exact  one. 

We  consider  condition  B  as  an  optional  one  because,  in  our  experience,  it  is  often  both 
hard  and  expensive  to  obtain  proper  experimental  data.  Furthermore,  these  data  might  be 
suitable  only  for  one  version  of  that  numerical  method  and  not  suitable  for  other  versions. 
Definition  6.1  is  worthy  of  some  discussion,  which  is  done  below  in  this  subsection. 

The  single  reason  why  we  have  introduced  Definition  6.1  is  a  substantial  challenge  in 
the  goal  of  constructing  such  a  numerical  method  for  a  CIP  which  would  provide  a  good 
approximation  for  the  exact  solution  without  an  advanced  knowledge  of  a  small  neighborhood 
of  this  solution  (section  1).  Because  of  this  challenge,  it  is  unlikely  that  the  desired  good 
approximation  for  the  exact  solution  would  be  obtained  without  some  approximations.  Thus, 
we  use  the  approximate  mathematical  model  M. 

The  main  requirement  of  Definition  6.1  is  that  this  numerical  method  should  provide  a 
sufficiently  good  approximation  for  the  exact  solution  x*,  regardless  on  any  a  priori  knowledge 
of  a  sufficiently  small  neighborhood  of  x*.  Furthermore,  it  is  important  that  one  should  have 
a  rigorous  guarantee  of  the  latter,  again  within  the  framework  of  the  model  M.  LInlike 
the  classical  convergence,  this  definition  does  not  require  that  lim  xn  =  x*.  Furthermore, 

n— >oo 

the  total  number  of  iterations  can  be  finite  and  should  be  chosen  as  a  result  of  numerical 
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studies.  We  remind  that  it  is  quite  often  in  the  field  of  Ill-Posed  Problems  that  the  number 
of  iterations  is  a  regularization  parameter,  see,  e.g.  pages  156  and  157  of  [19]. 

Therefore,  Definition  6.1  leaves  the  room  for  a  refinement  of  the  approximate  solution 
Xgiob  via  a  subsequent  application  of  a  locally  convergent  numerical  method.  In  other  words, 
the  room  is  left  for  a  two-stage  numerical  procedure  [5,  6,  7,  8].  In  accordance  with  the 
theory  of  Ill-Posed  Problems  [19,  34,  40],  some  a  priori  knowledge  about  the  exact  solution 
x*  should  still  be  in  place  of  course,  see,  e.g.  (3)  and  (84)-(87). 

Remark  6.1.  As  to  the  approximate  mathematical  model  M,  here  is  a  good  analogy.  It 
is  well  known  that  the  Huygens-Fresncl  optics  is  not  yet  rigorously  derived  from  the  Maxwell 
equations.  We  now  cite  some  relevant  statements  from  section  8.1  of  the  classical  book  of 
Born  and  Wolf  [12],  First,  “Diffraction  problems  are  amongst  the  most  difficult  ones  en¬ 
countered  in  optics.  Solutions  which,  in  some  sense,  can  be  regarded  as  rigorous  are  very 
rare  in  diffraction  theory.”  Next,  “because  of  mathematical  difficulties,  approximate  models 
must  be  used  in  most  cases  of  practical  interest.  Of  these  the  theory  of  Huygens  and  Fresnel 
is  by  far  the  most  powerful  and  is  adequate  for  the  treatment  of  the  majority  of  problems 
encountered  in  instrumental  optics.”  It  is  well  known  that  the  entire  optical  industry  nowa¬ 
days  is  based  on  the  Huygens-Fresnel  theory.  Analoguosly,  although  our  method  of  both 
the  current  and  the  above  cited  previous  publications  works  only  with  approximate  models, 
its  accurate  numerical  performance  has  been  consistently  demonstrated  in  the  above  cited 
publications. 

6.3  Our  approximate  mathematical  model 

Assuming  that  conditions  of  Lemma  3.1  hold,  Assumptions  1-3  below  mean  that  we  take 
into  account  only  the  first  term  of  the  asymptotic  behavior  of  the  function  s-2  In  w*  (x,s) 
at  s  — >  cx)  and  ignore  the  rest.  Assumption  4  means  that  we  assume  that  the  function 
q*(x,s)  satisfies  the  third  boundary  condition  (42).  Recall  that  the  latter  condition  is  an 
approximate  one. 

The  equation  for  the  function  q*  is  (see  formula  (85)  in  [4]) 


s 

A q*  —  2 s2Vg*  j  Vg*  ( x ,  r)  dr  +  2s 

1 

"<P  ' 

* 

yy 

a. 

-s 

S 

Is  J 

+  2s2Vg*  W* 


-  2sVR* 


Vg*  ( x ,  r)  dr  +  2s  (W*)2  =  0,  ( x ,  s)  E  x  (s,  s) . 


S 


(96) 


Our  approximate  mathematical  model  M  consists  of  the  following  four  assumptions: 
Assumptions.  We  assume  below  that: 

1.  The  exact  solution  of  our  Inverse  Problem  c*  ( x )  exists  and  satisfies  condition  (84). 
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2.  There  exists  a  function  p*  (. x )  E  H5  (12)  such  that  the  exact  tail  function  V*  ( x )  has 
the  form  (see  (90)) 

V*  (z,s)  V*  (x)  -  -  LiLjlA.  (97) 

s  s 

3.  For  s  >  s  >  s°  (d)  >  0  and  for  the  tail  function  V*  ( x )  in  the  form  (97)  the  function 
q*  (x,s)  ,(x,s)  E  12  x  [s,  s]  satisfies  condition  (85),  equation  (96).  In  addition  for  s  :=  s  the 
function  q*  (x,  s)  has  the  form 

q*  {x,s)  =  -^r--  (98) 

4.  For  s  E  [s,  s]  the  function  q*  (x,s)  satisfies  boundary  conditions  (88)  at  Ti  as  well  as 
the  third  boundary  condition  (42)  at  T2 

dnq*  (x,  s )  |r2=  s-2,  Vs  E  [s,  s] .  (99) 

In  addition  for  s  :=  s  the  function  q*  (x,s)  satisfies  the  third  boundary  condition  (42)  at  T3 


dsq*(x,s)  r3=  s  2. 


(100) 


It  follows  from  (99)  that 

dnq*n  |r2=  (snsn_i)_1 .  (101) 

Set  in  (96)  s  =  s.  Then,  using  (88)  and  (97)-(100),  we  obtain  the  following  approximate 
PDE  and  boundary  conditions  for  the  function  p*  (. x ) 

A p*  =  0  in  12,  p*  E  H5  (12) ,  (102) 

P*|ri  =  o  (®,s) ,  dnp*  |n  =  - s {x,s) ,  dnp*  |r2ur3=  -1-  (103) 


Boundary  conditions  (103)  are  over-determined  ones.  The  existence  of  the  solution  of  the 
problem  (102),  (103)  is  actually  assumed  because  conditions  (102),  (103)  are  derived  from 
Assumptions  1-4.  Let  functions  be  the  boundary  conditions  in  (42).  Suppose  that  for 

each  a  E  (0, 1)  there  exists  the  QRM  solution  p  =  p(x;  a)  of  the  following  boundary  value 
problem 

A p  =  0  in  12,  p  (x)  E  H5  (12) ,  (104) 

p|n  =  -sVo  (x,s) ,  dnp\Tl  =  -s2^ i  (x,s) ,  dnp  |r2ur3=  -1,  (105) 


see  Lemma  7.1  for  the  existence  and  uniqueness  of  the  function  p.  Then  we  choose  an  appro¬ 
priate  a  E  (0, 1)  and  in  the  iterative  process  of  subsection  4.3  we  set  the  first  approximation 
for  the  tail  function  as 

Vi,i  (x)  :=  V1}1  (x;  a)  :=  p^a\  (106) 

s 

To  show  that  equation  (102)  is  a  reasonable  approximation,  consider,  for  example,  the 
simplest  case  c  (x)  =  1.  Then  by  (14)  and  (90)  p*  (x)  =  —  \x  —  xo|  •  Hence, 

A p*  =  —  \x  —  xop1  ~  0  in  12,  (107) 
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as  long  as  the  source  xo  is  located  far  from  the  domain  hi,  as  it  is  the  case  in  many  applications. 
In  the  more  general  case  of  the  asymptotics  (36)  p*  (. x )  =  l  (x,xo)  >  \x  —  xq  \  .  Hence,  by 
an  analogy  with  (107),  assume  that  |AZ(x,  x0)|  <  CD1  (x,  x0) ,  x  G  fi,  where  C  >  0  is 
a  certain  constant  independent  on  x,x0.  Since  l(x,x0)  >  \x  —  x0|  ,  then  \Al(x,x0)\  < 
C  \x  —  Xop1 ,  x  G  Q.  Hence,  if  the  source  Xq  is  located  far  from  D,  then  similarly  with  (107) 
A l  (x,Xq)  ~  0  in  O.  The  function  l  (x,Xo)  satisfies  the  eikonal  equation  [38]  \Wxl  (x,  2:0)  |2  = 
c  (x) .  However,  one  cannot  consider  this  equation  in  our  case,  since  its  right  hand  side  is 
unknown.  There  is  of  course  a  well  known  the  so-called  “inverse  kinematic  problem”  which 
aims  to  recover  the  function  c  (x)  from  the  eikonal  equation  assuming  that  the  function 
l(x,x o)  is  known  for  all  x,Xq  G  dTl  [38].  However,  we  have  only  one  source  position  here. 
Consider  now  the  case  when  the  incident  plane  wave  is  originated  at  the  plane  {2:3  =  £30}  . 
The  function  w  corresponding  to  the  case  c  (x)  =  1  is  w  =  exp  (— s  (2:3  —  2:3,0!)  /(2s).  Hence, 
if  Q  C  {2:3  >  2:3,0}  and  c  (x)  =  1,  then  A p*  (x)  =  —A  (2:3  —  2:3,0)  =  0  in  H,  which  is  the  same 
as  (102). 

We  now  establish  uniqueness  within  the  framework  of  our  approximate  mathematical 
model.  Although  it  can  be  proven  under  less  restrictive  assumptions  imposed  on  functions 
q*,p*  than  ones  above,  we  are  not  doing  this  here  for  brevity. 

Lemma  6.1.  Suppose  that  above  assumptions  1-4  hold.  Then  there  exists  at  most  one 
function  q*  (x,  s )  for  (x,  s)  G  H  x  [s,  s] .  In  addition,  if  assuming  the  continuous  analog  of 
(93)  (as  in  [4,  5])  c*  (x)  =  Av*  (x)  +  s2  |Vv*  (2:)|2 ,  (x,  s)  G  H  x  [s,  s] ,  where  the  function  v 
is  the  same  as  in  (39)  with  (q,  V )  :=  (q*,  V *) ,  then  there  exists  at  most  one  function  c*  (x) . 

Proof.  We  outline  the  proof  only  briefly  because  it  is  simple.  Uniqueness  of  the  problem 
(102),  (103)  is  obvious.  Next,  substituting  in  (96)  the  function  V*  (x)  from  (97)  and  applying 
the  Carleman  estimate  of  Lemma  5.1  to  (96)  we  quickly  obtain  uniqueness  of  the  function 
q*(x,s).  The  s— integrals  are  not  a  problem  in  this  case  since  the  Carleman  estimate  is 
independent  on  lower  order  derivatives  and  it  can  be  integrated  with  respect  to  s  G  (s,  s) .  □ 

6.4  Estimates  for  the  tail  function 

Below  in  sections  6,7  B  =  B(s,d,x0 )  >  2  denotes  different  constants  depending  on  listed 
parameters.  We  do  not  indicate  its  dependence  on  the  domain  H,  because  is  as  in  (57)  in 
these  sections. 

Theorem  6.1.  Let  the  domain  H  be  as  in  (57)  and  the  source  xq  ^  fh  Let  the  function 
c*  (x)  satisfying  (84)  be  the  exact  solution  of  our  inverse  problem  as  defined  in  subsection 
6.1.  Let  0  <  s(d)  <  s  <  s  and  V*  (x)  be  the  exact  tail  function  as  in  (90).  For  each  function 
c  (x)  satisfying  condition  (9)  and  for  s  G  [s,  s]  let  w  (x,  s;  c)  :=  w  (x,  s )  be  the  solution  of  the 
problem  (11),  (12)  (Theorem  3.1).  Denote  V  (2 :;  c)  :=  V  (x)  :=  s~2  In  w  (x,  s;  c) .  Then  with 
a  constant  B  —  B  (s,  d,  xo) 

livr||c(n),||wic(n)  <  B,  (108) 

l|VV~VH|Mn)  +  ||AV-AV*|Un)  <  B|k-c*||Mn),  (109) 
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Proof.  Recall  that  x  =  (. x,y,z ) .  For  brevity  we  estimate  only  ||I4||c/j=;\  •  Estimates  of 
two  other  Erst  derivatives  are  similar.  By  (16)  and  (90) 

VO  VO  * 

\Vx\  =  —(x,s)  <  B\wx(x,s)\,  \v;\  =  -f(x,s)  <B\w*x(x,s) |.  (110) 

w  w* 

Theorem  3.1,  (14)  and  (30)  imply  that  for  £  =  (£i,  £2,  £3) ,  x  G  Cl,  b  (x)  —  c{x)  —  1 
wx  (x,  s)  =  w0x  (x,  s)  + 

s2  f\(x-^  ,  x  —  £,  \  J  /  -i  -i  -i  \ 


k-£iJ  k-£r 


exp  (— s  \x  —  £|)  b  (£)  w  (£,  s)  d£.  (Ill) 


Since  a;0  ^  fi,  then  functions  w0,w0x  do  not  have  a  singularity  for  x  E  Cl.  Hence,  (16)  and 
(111)  imply  that 


\wx  (x,  s)|  <  B  +  B  /  (s 


\x~0  \x  —  £|s 


exp  (— s  |x—  £|)  d£  <  B,  x  E  Cl.  (112) 


Hence,  (108)  follows  from  (110)  and  (112).  Next,  for  w  :=  w  —  w* 


vx  -  v:  = 


wx  vj; 

- -W  )  (z,s),  x  E  a 

W  ww*  / 


Hence,  by  (16)  and  (112) 

ll^^  — IIl2(T2)  —  B  ^||Vw||L2(n)  +  IMIl2(T2))  5 

B  (||Vw||  L2(  r3)  +  HHIl2(R3)) 

Let  c  =  c  —  c*.  Since  cw  —  c*w*  =  cw  +  cw *,  we  obtain  from  (11) 

Aw  —  s2c  (x)  w  =  s2cw*,  x  E  M3. 


(113) 


(114) 


Let  the  number  R  >  0  be  so  large  that  fl  C  Br  —  {|x|  <  i?}  .  Multiply  both  sides  of  (114) 
by  (—w)  and  integrate  over  Br.  Since  c(x)  =  0  outside  of  H,  we  obtain 

/  (|Vih|2  +  s2cw2)  (x,  s)  dx  —  f  fw^-\(x,s)dS  = 


— s2  J  (cw*w)(x,s)dx.(  115) 
n 

It  follows  from  (30)  and  (111)  that  Vw(x,s)  ,w(x,s)  E  L2  (M3)  and  the  second  term  in  the 
left  hand  side  of  (115)  tends  to  zero  as  R  — >  00.  Hence,  setting  in  (115)  R  — >  00,  we  obtain 


(|  Vw|2  +  s2cw 2)  (a;,  s)  dx  =  —s2  /  ( cw*w )  (x,  s )  dx. 


27 


Since  c  >  1,  then  s2cw2  (x,  s)  >  s2w2  (x,  s) .  Hence,  using  (16)  and  the  Cauchy-Schwarz 
inequality,  we  obtain 


w  (xi  s )  II  //qR3)  —  B  || 

cIIl2(0)  • 


(116) 


Next, 


AH -AH*  = 


Aw  V(u>  +  w*)  _  (Aw*  (Vw*Y(w  +  w*) 

- — - -Vw  - - - — — o - -  1  w 


w 


wr 


ww* 


(iuw* 


(x,s)(  117) 


Since  Aw*  (x,s)  =  s2  ( c*w *)  (x,s)  for  x  G  fl,  then  (16),  (112)  and  (117)  imply  that 


| AH  —  AH*|  <  B(\Aw\  +  |Vw|  +  \w\),x  G  fl.  (118) 

By  (16)  and  (114)  ||Aw||La(R3)  <  B  (||tZ?||i2(R3)  +  ||^Il2(q))  •  Hence,  (113),  (116)  and  (118) 
imply  (109).  □ 


7  Approximate  Global  Convergence  Theorem 

Assume  that 

s>l,  fih>l.  (119) 

Then  [4] 

max  {|Ai  n|  +  |A2  n|}  <  8s2.  (120) 

1  <n<N 

In  general,  embedding  theorems  are  valid  for  domains  with  sufficiently  smooth  boundaries. 
It  follows  from  Lemma  1  of  §4  of  Chapter  3  of  the  book  [35]  that  if  Q  is  a  rectangular 
prism,  then  any  function  /  G  Hk  ( Q )  can  be  extended  in  a  bigger  rectangular  prism  Qi 
D  Q,dQ  fl  dQi  =  0  as  the  function  f\  e  Hk  (Q\) ,  f\  (x)  =  f  (x)  in  Q  and  \\fi\\Hk^Q^  < 
Z  ||/||^fe(Q) ,  where  the  constant  Z  =  Z  (Q,Q i)  >  0.  Hence,  embedding  theorems  are  valid 
for  rectangular  prisms.  Hence, 

\\f\\c3(u)  <C\\f\\H5m,VfeH5(Q).  (121) 

Let  the  domain  fl  be  the  same  as  in  section  5.  Recall  that  fl*,  C  fl  for  x  e  (1/3, 1)  and 
fli  =  fl.  Following  the  construction  in  the  end  of  subsection  4.1,  we  assume  that 

P0  :=  x0  G  (1/3, 1) ,  c(x)  —  1  for  x  G  M3\flJ<D.  (122) 

Since  cn^,  (x)  Y  cn,k  (x)  for  x  G  fl^Xfl',  then  the  number  meas  (fl^0\fl')  can  be  considered 
as  a  part  of  the  error  in  the  data.  Hence,  we  assume  that  fl'  is  such  that 


meas  (fl^0\fl')  <  e/2, 


(123) 
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where  e  G  (0, 1)  is  sufficiently  small.  Since  by  construction  cny  (x)  ,  c*  (. x )  G  [1,  d] ,  Vx  G  M3 
and  cn,fe  (x)  =  cny  (x) ,  Vx  G  O',  then  by  (48),  (122)  and  (123) 

||cn,fc  —  c  <  llcn,fc  —  c  1 1 X/2  (s^')  +  —  llcn,fc  —  c  1^2(0')  4"  de.  (124) 

As  it  is  always  the  case  in  the  convergence  analysis  of  ill-posed  problems,  we  need  to 
connect  the  regularization  parameter  a  of  the  QRM  in  (56)  with  various  approximation 
errors:  the  level  of  the  error  a  in  the  data  (Lemma  7.1),  the  grid  step  size  h  in  the  s— direction 
and  the  number  £  in  (123).  In  addition,  to  iteratively  “suppress”  the  large  parameter  a~bl 
in  (83)  (to  ensure  the  stability  of  our  process),  we  need  to  impose  a  smallness  assumption 
on  the  length  (3  =  s  —  s  =  Nh,  where  N  >  1  is  an  integer.  For  a  number  x  >  0  let  {a:}° 
denotes  such  an  integer  that  x  —  {a:}0  G  [0, 1) .  Thus,  we  impose  the  following  conditions 

cr,E  G  (0,-s/o],  (125) 

h  =  y/a,  (3  :=  /3(a)  =  \/a{f(a)}°  :=  y/aN,  (126) 

where  the  function  /  (a)  is  monotonically  decreasing  for  a  G  (0, 1), 

/  (a)  >  0  for  a  G  (0, 1) ,  lim  /  (a)  =  oo  and  lim  ^  ^  =  0.  (127) 

a— >0+  a— >0+  In  (cc  x) 

Two  examples  of  the  function  /  (a)  are  fi  (a)  =  [In  (a-1)]'’ ,  r  =  const.  G  (0, 1)  and  fi  (a)  = 
In  (In  (a-1)) . 

Let  (qn,ki  cn,k,  Vn>k)  be  the  triple  computed  on  a  certain  step  of  our  iterative  process  of 
subsection  4.2.  Denote 


Qn,k  Qn,k  Qni  ^ n,k  C  ,  Vn,k  L n,k  . 

Similarly  for  qn,cn,  Vn.  Even  though  we  have  assumed  (for  brevity  only)  that  there  is  no  error 
in  functions  of  (94),  Lemma  7.1  and  Theorem  7.1  “allow”  error  to  be  present  in  functions 
-00  (x,s)  ,ij}*  (x,s)  in  (103),  see  the  first  two  sentences  of  the  proof  of  Lemma  7.1. 

Lemma  7.1  (estimate  of  V\^).  Let  the  domain  f 1  be  as  in  (57)  and  the  source  x0  ^  fb 
Let  Assumptions  1-4  of  subsection  6.3  hold  as  well  as  (125).  Let  T*  G  H5  (O)  be  a  function 
satisfying  boundary  conditions  (103).  Suppose  that  there  exists  a  function  T  G  H5  (O) 
satisfying  boundary  conditions  (105).  Let  the  number  a  G  (0, 1)  be  the  level  of  the  error  in 
the  function  T*  when  it  is  replaced  with  the  function  T  and 


ll*,-**llH»(n)  <  (128) 

l|p*llH»,n,  <  B.  (129) 

Let  the  function  p  —  p  (x;  a)  G  H5  (D)  be  the  unique  QRM  solution  of  the  problem  (104), 
(105)  which  is  guaranteed  by  Lemma  5.2,  and  the  tail  function  VL]  (x)  :=  V^i  (x;a)  has  the 
form  (106).  Then  for  every  a  G  (0, 1) 


l2(R) 


A  Vi 


l2<Si) 


<  By/a, 


HWullcfn) 


<  B. 


(130) 

(131) 
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Proof.  Note  that  the  trace  theorem,  (128),  (103)  and  (105)  imply  that 

114  (®,s)  -  4  (z,s)|4i(ri)  +  114  (x,s)  -  4  (^,s)IL2(ri)  <  Ca/s2. 

This  means  that  the  error  is  introduced  in  the  boundary  data  4  (x,  s) ,  4  (x,  s)  and  its  level 
is  proportional  to  a  G  (0,  ■\/a\  .  For  brevity  we  do  not  put  in  this  proof  the  dependence  of 
the  function  p  on  a.  The  existence  of  the  function  T*  follows  from  Assumptions  1-4.  Denote 
p(x)  =  (p  —  T )  (x)  —  (p*  —  T*)  ( x ) .  Then 

(A p,  Av)  +  a  [p,  v]  =  (AT*  —  AT,  An)  +  a  [T*  —  T,  v]  +  a  \p*,  v] , 

for  all  functions  v  G  H5  (D)  satisfying  zero  boundary  conditions  (105).  Setting  here  v  :=  p, 
and  using  (125),  (128),  and  (129),  we  obtain 

II^Plli2(n)  +  a  ll^llrr5(n)  —  a ^  ■  (132) 

Estimate  ||Ap||^2^  in  (132)  from  the  below.  We  have 

(Ap)2  =  (plx  +  ply  +  piz)  +  2  pxxpyy  +  2  pxxpzz  +  2  PyyPZZ,  (133) 

2  PxxPyy  9X  (^pxpyy  )  2 PxPyyx  (2pxPyy)  T  9y  (  2  PxPxy)  T  2  PXyi 

2 PxxPzz  i^PxPzz)  2 PxPzzx  ( — P :rP zz )  T  9Z  (  2 pxpXz )  T  2 Pxz 

and  similarly  for  2 pyypzz.  Integrate  (133)  over  fl  using  these  formulas  for  products.  Since  by 
(103)  and  (105)  dnp  |^=  0,  then  boundary  integrals  will  be  equal  to  zero.  Next,  use 

X 

Px  {x,  y,z)=  I  pxx  (£,  y,  Z )  df 

-1/4 

and  similar  formulas  for  py,pz.  Using  (132),  we  obtain 

cB 2  >  l|Apfla(n,  >  £  l|B”p||2Mn)  >C||Vp||^(n). 

M=2 

This,  (97),  (106),  (128)  and  (129)  imply  (130).  Next,  by  (121),  (125),  (128)  (129)  and  (132) 
l|Vp|lc(o)  —  C1  ||p||^5(n)  <  B.  The  latter  and  (106)  imply  (131).  □ 

Theorem  7.1  (approximate  global  convergence  property  of  our  algorithm).  Let  con¬ 
ditions  of  Lemma  7.1  hold  and  s  >  1.  Also,  assume  that  the  following  conditions  hold: 
(87)-(94),  (H9),  (122), (123)  and  (125)-(129).  Let  the  number  p  G  (x0, 1) ,  m  be  the  num¬ 
ber  of  inner  iterations  for  functions  qn,k,k  G  [l,m]  (section  f)  and  f  be  the  function 
in  (126),  (127).  Then  there  exists  a  constant  D  =  D  (s,d,x0,C*,  x0,  p)  >  1,  numbers 


30 


b\  =  b\  (s,d,xo,C*,  x0,  p)  G  (0,1/2),  62  =  1/2—  &i  defined  in  Theorem  5.1  and  a  sufficiently 
small  ay  =  ao  (s,  d,  xq,  C*,  x0,  p,  /,  m)  G  (0, 1)  sttc/i  that  the  following  estimates  are  valid 

^  ^  <  of2'2,  V  (n,  a)  G  [1,  TV]  x  (0,  a0) .  (134) 

II II  Z,2(0') 

Thus,  the  iterative  process  of  section  f  is  approximately  globally  convergent  of  the  level  ab 2/2 
in  the  norm  of  the  space  L2  (O')  (Definition  6.1). 

Proof.  Theorem  3.1  and  (48)-(51)  guarantee  the  existence  and  uniqueness  of  tails 
VU)k.  Note  that  because  of  (85)  and  (86),  the  estimate  (87)  is  not  changing  when  the 
number  TV  of  subintervals  of  the  interval  [s,  s]  increases  with  the  decrease  of  the  param¬ 
eter  a.  Let  (n,k)  G  [1 ,  TV]  x  [1  ,m\.  Assuming  that  the  desired  constant  D  is  found,  we 
hrst  estimate  D2nm+Aof2 .  Using  (126)  and  (127),  we  obtain  for  a  sufficiently  small  number 
a0  —  ao  (s,  d,  Xq,  C*,  x0,  p,  /,  m)  G  (0, 1)  and  for  all  a  G  (0,  a0) 

jj2nm+4 ap2  jyiNm+4^2  <- 

Dlexp  {—  In  (a-1)  [b2  —  (2mlnL>)  (/  (a) /In  (cU1))]  }  <  of 2^2.  (135) 

Below  in  this  proof  a  G  (0,  ay) .  Since  by  (48)  cnik{x)  >  1  in  then  ||cn,fc||L2(Q/)  > 
meas  (IT) .  Hence,  (135)  implies  that  it  is  sufficient  to  prove  that 


c„  -  c*||Mn,,  <  D2nmak,  V  (n,  a)  e  [1,  N }  x  (0,  a„) . 


(136) 


By  (43),  (91),  (94)  and  (101)  the  function  qn tk  is  the  QRM  solution  of  the  following 
problem 


A qn,k  -  Ai,n  (  X2  (x)  h.^Wqj  -  VVrhk  )  VqU)k  =  H,hk, 


Qn,k  | r i  b) zqn  k  |  r x  dnqn,k  | r 2  9) 


(137) 

(138) 
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Let  Qn  ( x )  be  the  last  line  of  (139).  We  now  estimate  this  function  using  (87),  (119),  (120), 
(126),  (127)  and  (131) 

||Qn||i2(n)  <  8 s2CC*V^f  («)  (C*  +  V^f  (a)  +  B)  <  ab\  n  G  [1,  N } .  (140) 

First,  we  estimate  q\^.  Denote  :=  AqiA  +  VVi,i ^qn,k-  We  obtain  from  (137)- 

(139)  that  the  function  qi  i  satisfies  boundary  conditions  (138)  as  well  as  the  following  integral 
identity  for  all  functions  v  G  H 5  (fl)  satisfying  (138) 

(Gu,i<fi,i,  Gitiv)  +  a  [gu,  v]  =  G'gi'u)  -  a  [<?*,  v] , 

Hhl  :=  -  A2,nVV1,1  (VLy  +  VD)  +  Qx.  (141) 

By  (108),  (119)  and  (120)  H-A^iVVi.illcm)  <  8 Bs2.  Hence,  using  Lemma  5.3,  Theorem  5.1 
and  (87),  we  obtain 


ll9i,i|lH5(n)  <  1/2 

1111-1  II H2(nXQ)  <  D  ( a  1 


#1,1 

Hhi 


L2(Q) 


+ 1  > 


l2(P) 


-I-  ( y. 
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Estimate  now 


#1,1 

#1,1 


l2<Si) 


.  By  (87),  (108),  (120),  (130),  (131),  (140)  and  (141) 


L2(Ci) 


<  8 s2C*Bab2  +  16 s2Bab2  +  ab 2  <  8 s2B  ( C *  +  3)  a 
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Assuming  that 


we  obtain 


#1,1 


(142) 

(143) 


L2(Q) 


D  >  8s2  B  (C*  +  6) , 

<  Dab 2.  Hence,  using  (142)-(144)  and  b2  =  1/2  —  b\,  we  obtain 


(144) 


Iki,i||tf5(n)  <  D2  {a  61  +  l)  ,  (145) 

1151,111^(0,)  <  O2  (a,/2-2t'  +  ak)  .  (146) 

Since  qlti  =  qXA  +  ql,  then  by  (87),  (121)  and  (145)  lead  to 

Il9i,illc.(n)  <D3(a-b'+2).  (147) 

We  now  estimate  HcyiH^,-^  .  It  follows  from  (46),  (47)  and  (92)-(94)  that 

ci,!  =  [-hAq1A  +  Ab,))  +  si  (-/iVgi.i  +  VHU)  [- KV  ( qn  +  ql)  +  V  (Vn  +  V*)] .  (148) 


By  (127)  and  (135) 


D2Nm+4ab2  <n~1. 


(149) 
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Hence,  (87),  (121),  (126),  (144)-(147)  and  (149)  imply  that 

'!llA9i,illL3(n,„).'>HV9i,i|li2(a,o)  <  2  D2a^<D3a^<a^N-',  (150) 

^I|V«iJc(j,)+2||V<,;||c(!5))  <  D3  (a*  +  ia1'2)  <  N~\  (151) 

Next,  by  (131)  ||  V  (V7,  ,  +  H*)||c/^\  <  B.  Hence,  using  (144)  and  (151),  we  obtain 

s?  ||  —hV  (qiti  +  ql)  +  V  (Vi,i  +  H*)||c/^\  <  D.  Hence,  (130),  (144),  (148)  and  (150)  imply 
that 

IPull^n.)  <  IPi,illl2(n,0)  <  (B  +  W-1)  (B  +  l)afc  <  D2ab\  (152) 

Hence,  (109),  (124),  (125),  (144)  and  (152)  lead  to 

VH12  +  AH12  <D3ab2.  (153) 

:  L2(n)  :  l2(Q) 

We  have  obtained  estimates  (145)-(147),  (150)-(153)  starting  from  the  estimates  (130), 
(131)  for  functions  Vi,i,  W,i,  V*.  Hence,  continuing  this  process  m  times,  using  q\  =  gl  m,  C\  = 
Citm  and  keeping  in  mind  that  by  (51)  Wi  =  Vi,m+i,  we  obtain  similarly  with  (150)-(153) 


'!llV9illl2(n,11),fc||A9.llt2(n,0)  <  (154) 

fe(l|V9illc(n)+2||V?;||c(!!))  <  N-\  (155) 

llCl_  c  II  L2(f2')  —  HCl'mllL2(n')  A  D  a2,  (156) 

VV2 1  +  AH21  <  D2m+1ab2.  (157) 

L2(Q)  :  L2(Q) 

Note  that  (156)  is  (136)  for  n  =  1.  LetG  [2,  At) .  Because  of  (154)-(157),  assume  that 

n— 1  n— 1  1 

^Ell^lkM^EllA^lkK)  <  (158) 

3  = 0  3=0 

n-1  /  \  _  1 

'>WllvyC(n)WIWIc(n))  £  V-  (159) 

3=0  '  ' 

VK,i  ,  N+  AW,  .  .  <  D2^m+1ab2  :=  IVi«\  (160) 

L/2yQ?cQj  L'2yQ?cQj 

\\cn-l  ~c  IIl2(0')  A  IlCn-l.mll^^,,)  A  D  ^  •*  CK  2.  (161) 


We  are  going  to  prove  now  (158)-(161)  and  (136)  for  n  :=  n  +  1.  Because  of  (137),  denote 

ii-i  \ 

y2  (x)  h  E  _  VW.i  )  v Qn,l ■ 

3=0  J 


^Qn,l  -^4l  ,n 


(162) 
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By  (94)  the  function  qn_\  satisfies  boundary  conditions  (138)  as  well  as  the  following  integral 
identity  for  all  functions  v  E  H5  (Q)  satisfying  boundary  conditions  (138) 

(Gn,iqn,i,  Gn> iv)  +  a  [qn> i,  v]  =  (hh,i,  GnAv^  -  a  [q*,  v] .  (163) 

Estimate  the  coefficient  at  V<fn, i  hi  (162).  Using  (108),  (120)  and  (159),  we  obtain 

(n—  1 

y2  (*)  h  ^  -  W„,i 

3=0 

In  terms  of  Theorem  5.1  an  important  feature  of  (164)  is  that  this  estimate  is  independent 
on  n,  implying  that  the  same  constants  D,b\,b2  =  1/2  —  b\  can  be  used  in  (165),  (166)  for 
all  n  E  [2,  N],  provided  that  (158)-(161)  are  valid  for  n  :=  n  +  1.  Thus,  using  Lemma  5.3, 
Theorem  5.1,  (144),  (163)  and  (164),  we  obtain 


<  16  Bs 


i-2 


(164) 


Ik™, ill H5(fi)  <  D  ^ 
?n,i||/f2(nxo)  <  D  fa  1 


Hn,  i 

Hn.  i 


l2(Q) 


+  1  > 


L2(S1) 


+  Of 
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(165) 

(166) 


Hence,  using  (85),  (120),  (121),  (139),  (140),  (144),  (158)-(161),  Theorem  6.1  and  that  B  >  2, 
we  obtain 


Hn.  1 


Mfi) 


<  8sz 


n  —  1 


-jy—  ck2  +  Dn_iab2  )  C*  +  8s2 — — 


n-  1  b2  (n  -  1 


-a 


N 


+  2  B 


-\~8s2  Dn_\oi 


b2 


2  (n  -  1) 
N 


+  2  B)+a 


,b2 


<  8s2Dn_1ab2  (3  +  4 B  +  3C*/2)  <  DDn_xab2. 
Hence,  (121),  (149),  (160),  (165)  and  (166)  imply  that 


h 


^11  ^9n,i||(7(n)  h"  2  ||Vgnil 


<  D3Dn_1ab2  <  D2Nm+4ab2  <N~\ 


(167) 


h  \\qn,i\\m^0)  <  D  (DDn- ia2b2  +  ab2+1/2)  <  D^+V62  <  ab2N~\  (168) 

We  obtain  similarly  with  (148) 

n- 1  \  /  n- 1 

~hAqn<  1  —  h  A qj  +  Al/t,i  )  +  |  —h\7qnj  1  —  h  V^-  +  VV/p 

3=0  )  \  3=0 

n—1 

■  -hV  (qn,  1  +  o  -  h  v  +  ?;)  +  V  (K,1  +  V*) 

3=0 
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Hence,  using  (144),  (158)-(160),  (167)  and  (168),  we  obtain 

IIc™.iIIl2(^0)  -  (jv^2  +  Dn-l(yb 2)  1  +  (iV  +  jB)  “  DDn-iah'2. 

Hence,  (109),  (124)  and  (125)  imply  that 

VK,2  +  AK,2  <D2Dn_xah\ 
l2(o)  l2(h) 

Similarly  for  k  —  1,  ...,m 

HCn4'llL2lo  )  —  D2k  1Dn_iah2,  S7Vnk+i  +  AVnk+i  <  D2k Dn_iab2 . 

L2(0)  L2(Q) 

Hence,  similarly  with  the  above,  we  obtain  that  estimates  (167),  (168)  are  valid  for  functions 
qn.k,  gnj,:.This  implies  the  validity  of  (158)  and  (159)  for  n  :=  n  +  1.  Similarly 

IK,*  -  c*IP=Coo  <  IK,*llla(n,„)  <  B2lA,-iate  =  fl2("-1)m+aak,fc  e  [l,m], 

,  +  AV',  „,  <  ,o'’'  =  D-""  iq'A 

L2(Q)  L2(n) 

The  last  two  estimates  establish  (160)  and  (161)  for  n  :=  n  +  1.  □ 

8  Numerical  Studies 

8.1  Main  discrepancies  between  the  convergence  analysis  and  the 
numerical  implementation 

8.1.1  Discussion  of  the  topic  of  discrepancies 

As  to  the  general  philosophy  of  numerical  solutions  of  challenging  nonlinear  problems,  both 
well-posed  and  ill-posed  ones,  it  is  well  known  that  some  discrepancies  between  the  conver¬ 
gence  analysis  and  numerical  implementations  are  almost  inevitable.  The  main  reason  is 
that  because  of  the  complicated  structure  of  those  problems,  the  theory  usually  can  grasp 
only  a  part  of  numerical  studies  rather  than  all  aspects.  In  general,  given  a  complicated 
nonlinear  ill-posed  problem,  it  is  almost  impossible  to  obtain  accurate  computational  results 
if  following  the  theory  exactly.  Thus,  it  is  well  known  that  computations  are  usually  carried 
out  under  less  restrictive  conditions  than  ones  imposed  by  the  convergence  analysis. 

For  example,  as  a  rule,  constants  in  convergence  theorems  are  significantly  over/  under¬ 
estimated  (maybe  with  the  only  exception  of  a  few  very  simple  linear  problems).  In  addition, 
the  theory  works  for  a  continuous  case  meaning  an  infinitely  dimensional  space.  On  the  other 
hand,  in  real  computations  one  always  works  with  a  discrete  case  of  the  FEM/FDM  in  a  finite 
dimensional  functional  space.  Furthermore,  in  the  case  of  an  ill-posed  problem  the  minimal 
spatial  grid  size  hsp  of  the  FEM/FDM  often  serves  as  an  implicit  regularization  parameter 
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for  the  discrete  case  (naturally  hsp  is  different  from  the  above  step  size  h  in  the  s— direction). 
Thus,  hsp  cannot  be  taken  too  small:  otherwise  that  functional  space  would  effectively  be¬ 
come  an  infinitely  dimensional  one.  For  example,  in  our  numerical  testing  below  we  have 
used  hsp  =  0.125  to  minimize  the  QRM  functional  (56)  via  the  FDM.  However,  our  attempt 
to  decrease  it  to  hsp  =  0.0625  has  significantly  worsened  results.  The  bound  from  below 
for  hsp  imposes  an  upper  bound  on  the  dimension  of  the  corresponding  finite  dimensional 
functional  space.  Thus,  already  over/under-estimated  constants  of  convergence  theorems  for 
the  continuous  case  become  even  more  over/under-estimated  for  the  finite  dimensional  case 
with  an  upper  bound  on  its  dimension. 

For  example,  Theorem  7.1  imposes  smallness  conditions  on  the  parameter  a  as  well  as 
on  other  parameters  in  (125)-(127),  which  is  an  “ideal  limiting  case”  for  our  problem.  We 
note  that  a  similar  situation  quite  often  occurs  when  one  proves  convergence  theorems  about 
nonlinear  problems,  both  well-posed  and  ill-posed  ones,  assuming  that  some  parameters  are 
“sufficiently  small” .  Despite  this  assumption,  it  is  often  the  case  in  real  computations  that 
those  parameters  are  not  too  small.  Another  example  is  an  accurate  imaging  from  experi¬ 
mental  data  with  a  large  noisy  component  in  [8,  27]  despite  the  fact  that  the  convergence 
analysis  of  [4,  5]  requires  a  sufficiently  small  error  in  the  input  data.  Figures  3-5  in  [27]  as 
well  as  Figure  2-a  in  [8]  give  an  idea  about  the  amount  of  the  noise  in  those  data. 

Having  stated  the  above,  we  nevertheless  point  out  that  the  convergence  analysis  is  ob¬ 
viously  a  very  important  element  of  numerical  studies.  This  analysis  provides  an  analytical 
guidance  for  computations,  qualitatively  explains  numerical  results  and  ensures  that  at  least 
in  an  “ideal  limiting  case”  the  convergence  is  guaranteed.  In  addition,  regardless  on  some 
deviations,  numerical  implementations  are  usually  guided  by  convergence  theorems.  Thus, 
a  certain  reasonable  balance  between  the  convergence  analysis  and  the  numerical  implemen¬ 
tation  should  take  place. 

In  particular,  speaking  about  our  problem,  we  believe  that  various  conditions  in  our 
convergence  analysis  might  be  relaxed  if  working  with  a  discrete  case  with  a  lower  bound 
imposed  on  hsp,  hsp  >  ho  =  const  >  0.  To  analyze  the  QRM  in  this  case,  one  needs  to  apply 
a  discrete  version  of  the  Carleman  estimate  of  Lemma  5.1,  see,  e.g.  [13,  25]  for  discrete 
Carleman  estimates  for  elliptic  operators.  The  latter,  in  turn  should  likely  result  in  relaxed 
conditions  of  a  discrete  analog  of  Theorem  7.1.  This  will  be  one  of  topics  of  our  future 
research. 

8.1.2  Main  discrepancies 

The  first  main  discrepancy  is  with  regard  to  Lemma  3.1  about  an  easily  verifiable  sufficient 
condition  of  the  regularity  of  geodesic  lines.  In  general,  such  a  condition  is  unknown,  except 
of  the  trivial  case  when  the  function  c  ( x )  is  close  to  a  constant.  On  the  other  hand,  the 
authors  are  unaware  about  any  reasonable  results  for  CIPs  for  hyperbolic  PDEs  without 
either  the  assumption  of  the  regularity  of  geodesic  lines  or  a  somewhat  close  assumption. 
We  verify  the  asymptotic  behavior  (36)  computationally,  see  subsection  7.2  of  [4], 

The  second  main  discrepancy  is  that  we  replace  in  our  computations  a  |M|#5(n)  in  (56) 
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1 1 2 

with  a  ,  because  the  latter  is  simpler  to  implement  numerically.  One  of  the  rea¬ 

sons  why  this  works  computationally  might  be  because  of  the  above  arguments  about  finite 
dimensional  spaces,  since  all  norms  are  equivalent  in  such  spaces. 

The  third  main  discrepancy  is  that  we  conduct  computations  for  the  case  when  the  point 
source  in  (2)  is  replaced  with  the  plane  wave.  This  is  because  our  current  software  works  with 
the  plane  wave  only  and  it  would  take  a  time  consuming  effort  to  make  it  work  for  the  point 
source.  In  addition,  the  case  of  the  plane  wave  is  reasonable  for  our  target  application  to 
imaging  of  plastic  land  mines,  since  the  wave  radiated  by  a  point  source  effectively  becomes 
a  plane  wave  when  that  source  is  located  far  from  the  domain  of  interest.  We  have  chosen 
the  point  source  in  (2)  with  the  goal  to  use  Lemma  3.1  which  actually  can  be  derived  from 
the  construction  of  the  fundamental  solution  of  a  general  hyperbolic  equation  in  Theorem 

4.1  of  [38].  Other  than  this,  the  above  technique  can  be  extended  to  the  case  of  the  plane 
wave  and  this  might  be  one  of  topics  of  our  future  research. 

When  this  manuscript  was  ready  for  submission,  M.V.  Klibanov  and  N.  Pantong  (a 
coauthor  of  [30])  have  implemented  the  case  of  the  point  source.  Computational  results  were 
almost  the  same  as  ones  below.  However,  keeping  in  mind  considerations  of  brevity  as  well 
as  the  necessity  to  submit  a  separate  paper  with  Pantong  as  a  coauthor,  the  authors  have 
decided  not  to  present  the  point  source  case  here. 

The  fourth  main  discrepancy  is  that  we  conduct  computations  for  the  2-D  case  since 
this  case  is  both  easier  to  implement  numerically.  Another  important  factor  is  that  the  2-D 
case  is  faster  to  compute.  The  3-D  case  is  quite  feasible  [30].  Indeed,  the  3-D  image  in  [30] 
was  computed  from  the  forward  problem  data  simulated  in  3-D.  At  the  same  time,  twenty 
four  2-D  inverse  problems  in  corresponding  2-D  cross-sections  were  solved  simultaneously  on 
twenty  four  processors  and  the  3-D  image  was  formed  this  way. 

8.2  A  simplified  model  of  imaging  of  plastic  land  mines  and  some 
details  of  the  numerical  implementation 

We  outline  both  topics  of  this  subsection  only  briefly  here  referring  to  [30]  for  details.  In 
the  spirit  of  sub-subsection  8.1.1,  we  have  made  some  simplifications  in  our  mathematical 
model.  The  first  main  simplification  of  our  model  is  that  we  consider  the  2-D  instead  of  the 
3-D  case.  Second,  we  ignore  the  air/ground  interface,  assuming  that  the  governing  PDE  is 
valid  on  the  entire  2-D  plane.  Indeed,  our  experience  of  working  with  the  experimental  data 
of  [31]  has  indicated  that  the  reflection  from  the  air /ground  interface  can  be  removed  from 
the  data  via  a  new  data  pre-processing  procedure.  The  idea  of  a  similar  data  pre-processing 
procedure  can  be  found  in  [8,  27].  Let  the  ground  be  (x  —  (x,z)  :  z  >  0}  C  M2.  Suppose 
that  a  polarized  electric  held  is  generated  by  a  plane  wave,  which  is  initialized  at  the  line 
{z  =  z°<0,i6l}  at  the  moment  of  time  t  —  0.  As  it  was  pointed  out  in  section  2,  the 
following  hyperbolic  equation  can  be  derived  from  the  Maxwell  equations 


£r(x.)utt  =  A u,  (x,  t)  G  R2  x  (0,  oo) , 
u  (x,  0)  =  0,  ut  (x,  0)  =  5  (z  —  £°)  , 


(169) 

(170) 
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where  the  function  w(x,  t)  is  a  component  of  the  electric  field  and  £r  (x)  is  the  spatially 
distributed  dielectric  constant.  We  assume  that  the  function  er  (x)  satisfies  conditions  (3), 
(4)  (in  2-D).  Let  the  function  tc0  (z,  s )  =  exp  (— s  \z  —  z0\)  /(2s)  be  the  one  which  correspond 
to  the  Laplace  transform  (9)  of  the  incident  plane  wave  with  £r(x)  =  1.  Applying  (9)  to 
(169),  we  obtain  the  following  analog  of  the  problem  (11),  (12) 

A w  —  s2£r(x)w  =  —  <5  —  z°)  ,  s  >  s°  =  const.  >  0,  x  G  IR2,  (171) 

lim  ( w  —  w0)  (x,  s)  =  0.  (172) 

|x|— >oo 

It  is  well  known  that  the  maximal  depth  of  an  antipersonnel  land  mine  does  not  exceed 
about  10  cm=0.1  m.  So,  we  model  these  mines  as  small  squares  with  the  0.1  m  length  of 
sides,  and  their  centers  are  at  the  depth  of  0.1  m.  We  set 

12  =  (x  =  (x,  z )  G  (—0.30,  0.30)  m  x  (0,  0.6)m}  .  Introducing  dimensionless  spatial  vari¬ 
ables  x'  =  x/  (0.1m)  without  changing  notations,  we  obtain  fl  =  (—3,3)  x  (0,6) ,  sides  of 
small  squares  modeling  mines  equal  1  and  depths  of  their  centers  are  at  {z  =  1} .  We  took 
flp0  =  (—3,3)  x  (0,3)  (subsection  4.1).  Tables  of  dielectric  constants  [39]  show  that  er  =  5 
and  er  =  22  in  the  dry  sand  in  the  trinitrotoluene  (TNT)  respectively,  meaning  22/5  =  4.4  of 
the  mine/background  contrast.  Hence,  considering  new  parameters  £'r  =  ey/5,  s'  =  s-0.1  •  \/5 
and  not  changing  notations,  we  obtain 

£r  (dry  sand)  =  l,£r  (TNT)  =  4.4.  (173) 

Because  of  (173)  and  (3),  we  impose  £r  (x)  G  [l,8],er(x)  =  1  outside  of  the  rectan¬ 
gle  Qp0 .  We  have  modified  our  algorithm  of  section  4  via  considering  functions  v  ( x ,  s )  = 
s~2  In  (■ w/wq )  ,q(x,  s )  =  dsv  (ay  s )  instead  of  v  (ay  s )  =  s~2  In  ( w ) ,  q  (ay  s )  =  dsv  (ay  s ) .  This 
has  resulted  in  obvious  modifications  of  above  equations.  A  slight  modification  of  Theorem 
7.1  can  be  proven  for  this  case.  We  have  observed  in  our  computations  that  at  the  mea¬ 
surement  side  Tx  =  {x  G  (—3,  3) ,  z  =  0}  of  the  above  square  fl  the  ratio  ( w/w0 )  (x,  0,  s)  ~  1 
for  s  >  1.2  which  means  a  poor  sensitivity  of  Tx  to  the  presence  of  abnormalities  for  values 
of  the  pseudo  frequency  s  >  1.2.  The  best  sensitivity  was  for  s  G  (0.5, 1.2) .  Hence,  one 
should  expect  that  the  modified  tail  function  V  (x,  s)  =  V  (x,  s)  —  s~2  In  w  (x,  s)  ~  0  for 
s  >  1.2  at  least  for  x  close  to  Tx.  Hence,  we  have  chosen  s  =  1.2  and  s  G  [0.5, 1.2]  :=  [s,  s] . 
However,  if  we  would  work  in  the  original  domain  making  spatial  variables  dimensionless 
as  x"  =  x/  (lm)  ,  then  s"  =  yfb, s  implying  that  s"  =  12,  which  can  be  considered  as  a  large 
pseudo  frequency.  The  latter  shows  that  in  practical  computations  the  above  notion  of  suffi¬ 
ciently  large  s  is  actually  a  conditional  one  and  depends  on  particular  ranges  of  parameters 
at  hands  (in  fact,  this  is  in  the  spirit  of  subsection  8.1). 

To  simulate  the  data  at  {z  =  0}  for  the  inverse  problem,  we  have  solved  the  problem 
(171),  (172)  via  the  FDM  for  a  number  of  points  sn  G  [0.5, 1.2],  For  comparison,  we  have 
also  solved  once  the  problem  (169),  (170)  and  have  applied  the  Laplace  transform  (9)  then. 
Imaging  results  were  almost  the  same.  We  have  added  the  random  noise  of  the  5%  level 
to  the  simulated  Dirichlet  data  at  Tx  via  wa  (ay,  0,  sn )  =  w  (ay,  0,  sn)  (1  +  aun) ,  a  =  0.05, 
where  {ay}  are  grid  points  of  the  FDM  for  the  forward  problem  solution  and  u  G  (—1, 1) 
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is  a  random  variable.  To  calculate  the  derivative  ds  [s^2  In  ( wa/wo )  (xi,  0,  sn)]  (to  obtain  the 
boundary  data  for  q(xi,sn )),  we  have  smoothed  first  values  of  wa  (xj,0,  sn)  with  respect  to 
s  via  cubic  B— splines  similarly  with,  e.g.  [22],  Next,  we  have  used  finite  differences  to 
calculate  the  desired  derivative. 

The  modified  functional  (56)  was  written  in  the  FDM  form.  Its  minimization  was  per¬ 
formed  with  respect  to  values  of  the  function  u  at  grid  points  via  the  conjugate  gradient 
method.  We  have  chosen  a  =  0.04  and  the  spatial  grid  step  size  hsp  =  0.125.  First,  we  have 
solved  the  problem  (104),  (105)  via  the  QRM  and  thus  have  calculated  the  first  tail  Vjj  {x) 
in  (107).  Next,  we  have  continued  as  in  section  4  with  m  :=  10.  The  above  algorithm  has 
provided  us  with  the  function  £r.gi0b  (x)  (Definition  6.1). 

We  have  observed  in  our  computations  that  the  above  algorithm  can  accurately  image 
locations  of  mine-like  targets.  However,  values  of  the  function  £r,giob  (x)  near  points  of  local 
maxima  were  not  imaged  accurately.  Thus,  similarly  with  [5,  6,  7,  8,  30],  we  have  applied 
a  two-stage  numerical  procedure,  see  comments  in  subsection  1.1.  While  the  first  stage  was 
the  one  described  above,  on  the  second  stage  we  have  minimized  the  Tikhonov  regularization 
functional  via  the  gradient  method  taking  the  function  £r,giob  (x)  as  the  starting  point.  Since 
the  second  stage  is  not  a  focus  of  the  current  paper,  we  refer  to  [30]  for  a  detailed  description. 
Our  attempt  to  use  the  gradient  method  without  the  first  stage  did  not  lead  to  accurate 
results.  The  latter  indicates  the  importance  of  the  first  stage. 

8.3  Numerical  results 

In  this  subsection  we  present  two  numerical  examples  of  the  performance  of  our  numerical 
method  for  computationally  simulated  data  using  the  above  mathematical  model  of  imaging 
of  plastic  antipersonnel  land  mines. 

Test  1.  We  test  our  numerical  method  for  the  case  when  plastic  land  mines  are  modeled 
by  two  small  squares  with  the  same  size  sz  —  1  of  their  sides.  Also,  £r  =  4.4  in  both  of 
them  and  £r  —  1  everywhere  else,  see  (173)  and  Figure  8.1-a.  Centers  of  these  squares  are 
at  (x*,2i)  =  (—1.5, 1.0)  and  (xl,z%)  =  (1.5,1).  However,  we  do  not  assume  a  priori  in  our 
algorithm  neither  the  presence  of  these  squares  nor  a  knowledge  of  £r  (x)  at  any  point  of  the 
rectangle  fh  Figure  8.1-b  shows  the  image  of  £r^i0b  (x).  Locations  of  both  local  maxima  are 
imaged  well.  However,  values  of  those  maxima  are  2.8,  which  is  about  36%  off  the  correct 
one.  Figure  8.1-c  displays  the  final  image  after  applying  the  second  stage  of  our  two-stage 
numerical  procedure.  Locations  of  both  targets  are  imaged  well.  In  addition,  maximal  values 
of  the  computed  coefficient  £rjinai  in  both  targets  are  max  ( £r, final )  —  4.4,  which  is  accurate. 
Likewise,  £r, final  (x)  =  1  outside  of  those  spots,  which  is  also  accurate. 

Test  2.  We  test  our  numerical  method  for  the  case  of  the  same  two  squares  as  ones 
above.  Now,  however  we  have  £r  =  6  in  the  left  square,  £r  =  4  in  the  right  square  and  £r  —  1 
everywhere  else,  see  (173).  Again,  we  have  obtained  accurate  images,  see  Figure  8.2  and  the 
legend  for  it. 

Although  mine-like  targets  are  symmetrically  located  in  both  tests,  values  of  the  function 
£r  (x)  are  not  symmetric  in  Test  2.  A  similar  image  quality  for  the  case  of  asymmetrical 
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locations  of  two  targets  was  obtained  in  [30]. 


a) 


o  o 


b) 


o  o 


Figure  1:  Test  1.  The  plane  wave  falls  from  the  top  and  the  data  are  collected  on  the  top  boundary, 
a)  Correct  image  of  two  mine-like  targets,  er  =  4.4  in  both  small  squares  and  er  =  1  everywhere 
else,  see  (173).  b)  The  image  of  the  function  £r,giob  (x) .  While  locations  of  targets  are  imaged 
accurately,  values  of  local  maxima  are  2.8,  which  is  36%  off  the  correct  one.  c)  The  final  image 
of  the  function  £r, final  (x)  after  applying  the  second  stage  of  the  two-stage  numerical  procedure. 
Locations  of  both  mine-like  targets  are  accurately  imaged.  Likewise,  max  (£r, final)  =  4.4  in  both 
spots  and  £r, final  =  1  outside  of  these  spots,  which  is  also  accurate. 
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